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ABSTRACT 


The  BOXPLX  of  a  user  specified  cost  function  is  evaluated 
as  a  control  system  design  method.   The  emphasis  is  on  compen- 
sator design  for  multiple-input-multiple-output  plants,  but  the 
technique  should  be  applicable  to  any  control  system  with 
adjustable  parameters. 

The  engineer  must  define  the  plant  dynamics,  inputs  and 
outputs,  the  desired  output  responses  for  the  specified  inputs, 
weighting  factors  for  the  performance  measure,  the  type  of  com- 
pensation to  be  used,  the  free  parameters,  and  initial  values 
and  constraints  on  the  free  parameters . 

The  program  gives  a  direct  search  of  the  optimum  values, 
within  the  specified  constraints ,  for  the  free  parameters . 
Evaluation  of  the  results  is  given  in  graphical  plots  of  the 
output  time  responses  including  the  desired  output  response  as 
well  as  the  compensated  response. 

The  Transfer  Matrix  Method  is  used  to  provide  an  analytical 
check  on  the  accuracy  of  the  method  and  the  procedure  is  illus- 
trated with  a  two- input- two-output  system  example. 
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I.   INTRODUCTION 

As  control  systems  become  more  and  more  complex,  design 
engineers  are  faced  with  an  increasingly  difficult  problem  in 
the  development  of  compensators.   Single-input-single-output 
linear  systems  can  be  compensated  using  Bode  plots  or  root 
locus  plots  and  trial  and  error  methods  to  meet  the  system 
specifications  [1] .   When  compensating  multiple-input-multiple- 
output  linear  systems,  use  can  be  made  of  transfer  matrix 
techniques  described  by  Ogata  [2]  to  achieve  total  decoupling 
and  the  desired  output  responses.   Since  decoupling  during 
the  transient  period  of  real  systems  may  require  controls 
which  exceed  physical  constraints  imposed  by  the  system,  it  is 
sometimes  desirable  to  decouple  during  steady  state  conditions 
and  allow  coupling  during  the  transient  period  [3]  . 

Non-linear,  single-input,  single-ouput  systems,  wherein  the 
non-linearity  can  be  lumped  into  one  element,  can  be  analyzed 
and  compensated  effectively  using  the  phase  plane  method  (for 
second  order  systems)  and  describing  function  analysis  [2,4]  . 
Seme  of  these  ffiethods require  extensive  mathematical  manipulations 
and  become  rather  unwieldy,  but  do  produce  the  desired  results. 
In  ail  of  the  analysis  and  design  approaches  mentioned  above, 
any  one  of  several  well  established  simulation  programs  can  be 
used  to  evaluate  the  design. 

To  reduce  the  tedium  involved  in  the  design  of  compensators , 
frequency  domain  and  time  domain  computer  optimization  schemes 
have  been  developed,  which  will  set  free  parameters  in  the 
compensator  to  yield  a  system  response  which  most:  nearly 


approximates  a  reference  response  specified  by  the  user.   Lima 
used  this  approach  to  address  the  problem  of  controller  design 
for  ships  engaged  in  underway  replenishment  [5].   MacNamara 
applied  this  method  to  an  autopilot  design  [6  3,  and  Vines 
developed  a  generalized  program  for  compensator  optimization 
in  single-input-single-output  systems  [7]. 

The  following  chapters  present,  a  generalized  discussion  of 
transfer  matrix  analysis  and  compensation  of  multivariable 
systems,  optimization  techniques,  and  performance  measures. 
A  two-input-two-output  system  is  then  compensated  using  transfej 
matrix  techniques,  and  using  parameter  optimization.   The 
results  are  then  compared. 


II.   MULTIVARIABLS  SYSTEM  COMPENSATION 

A.   TRANSFER  MATRIX-ANALYTICAL  METHODS 
1.   Introduction 

A  linear,  multiple-input-multiple-output  or  multivariable 
system  can  be  described  by  a  transfer  matrix,  which  is  simply 
an  extension  of  the  transfer  function  concept.   Consider  the 
two-input-two-output  open  loop  plant: 


RMS] 


R2(S) 


cils] 


Figure  2-1   Multivariable  Plant 

The  Block  Diagram  of  the  plant  can  be  represented  in  the 
general  sense  as  shown  in  Figure  2-2. 


RMS) 


R2(S) 


CMS] 


C2(S) 


Figure  2-2   Generalized  Multivariable  Plant 
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Algebraic  equations  for  the  outputs  in  the  independent  variable 
S,  can  be  written  directly: 


Cl(s)  =  GplKs)Rl(s)  +  Gpl2(s)  R2(s) 
C2(s)  =  Gp21(s)Rl(s)  +  Gp22(s)  R2(s) 


(1) 
(2) 


If  the  inputs  Rj (s)  and  the  outputs  C. (s)  are  each  considered 
as  vectors  of  one  dimension  (i,  j  =  1,  2)  the  equations  can  be 

written  in  matrix  form: 


Cl(s) 
C2(s)j 


Gpll(s) 
_Gp21(s) 


Gpl2(s) 


Gp22(s)_      _R2(s)_ 


RKs) 


(3) 


The  number  of  inputs  and  outputs  can  be  increased,  and  need  not 
be  equal  so  that  in  the  general  case: 


Cl(s) 

Gpll(s) 

C2(s) 

Gp21(s) 

1 

= 

1 

1 

1 

1 

' 

1 

Ci(s) 

Gpil(s) 

Gpl2(s)  Gpij(s) 


Gpij (s) 


RKs) 
R2(s) 


(4) 


R.  (s) 


-J    LJ    _i 


or 


C(s)   =  G  (s)R(s) 


(5) 


where  Gp(s)  is  the  transfer  matrix  (ixj)  of  the  plant. 

2.   Feedback  Compensation 

Transfer  matrices  can  be  used  to  describe  the  compen- 
sation of  multivariable  systems  as  well  as  the  plant.   Consider 
the  two-input-two-ouput  plant  above  with  feedback  compensation 
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as  shown  in  Figure  2-3. 


RKS) 


R2(S)  \£ 


U1[S] 


H  ,2  (S] 


<r 


H21(S)<" 


U2(S) 


ci(s) 


> 


C2(S) 


> 


Figure  2-3 
Multivariable  Plant  with  Feedback  Compensation 

Algebraic  equations  for  R,(s)  and  R9(s)  can  be  written: 
Rl(s)  =  Ul(s)  _  H  11(s)  cl(s)  _  H12(s)  C2(s) 
R2(s)  =  U2(s)  -  H  21(s)  Cl(s)  -  H22 (s)  C2(s) 


(7) 


with  C,  R,  and  U  two  element,  one  dimension  vectors  these 
equations  can  be  rewritten: 


Rl(s 


,.n 


R2(s 


UKs) 


U2(s 


Hll(s)     H12(s) 


H21(s 


H22  (s 


CI  (3 


C2(s 


(8) 
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or  in  the  general  case: 

P.(s)  =  U(s)   -   H(s)C(s)  (9) 

where  H(s)  is  the  transfer  matrix  (jxi)  of  the  compensation. 
From  the  uncompensated  plant: 

C(s)   =   G  (s)R(s)  (10) 

Substituting  equation  (9)  for  R(s)  in  Equation  (10): 

C(s)   =   G  (s)   [  U(s)  -  H(s)  C(s)]  (11) 

-P       - 

Rearranging: 

-[I  +  G  (s)H(s)  ]C(s)  =  G  (s)U(s)  (12) 

Since  Gp(s)  is  an  ixj  matrix,  K(s)  is  a  jxi  matrix  the  product, 
Gp(s)H(s),  is  a  square,  ixi  matrix.   Assuming  that  [I  +  Gp(s)H(s)] 
is  non-singular  the  inverse  can  be  taken. 
Premultiplying  both  sides  of  equation  (12)  by  [I+G  (s)H(s)]   : 

C(S)  »  [I  +  G  (s)H(s)]"1G  (s)U(s)  (13) 

C(s)  =   G(s)U(s)  (14) 

The  compensated  system  transfer  matrix  ( ix j )  is  then: 

G(s)  =  [I  +  G^(s)H(s)]  "1Gr^(s)  (15) 

-P  -P 

The  plant  with  feedback  compensation  can  be  depicted  more 
clearly  as  shown  in  Figure  2-4. 
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CIS) 


^> 


Figure  2-4 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Feedback  Compensation 


To  determine  the  feedback  compensation,  H(s),  required  to 

achieve  design  specifications,  the  dynamics  of  the  plant  must 

be  known  and  put  into  transfer  matrix  form,  G  (s) .   The  desired 

closed  loop  dynamics  must  then  be  reflected  in  the  closed  loop 

transfer  matrix,  G(s) .   Then  through  matrix  manipulation  the 

compensation,  H(s) ,  can  be  computed  in  closed  form  provided 

that  G(s)  and  G  (s)  are  scuare  and  non-singular. 
-        -p 

Premultiolying  both  sides  of  equation  (15)  by  [I+G  (s)H(s)] 
*■-*.-'  ^  -iL„-p..j 

and  then  postmultiplying  by  G(s)~  : 


I  +  6  (s)H(s)  =  G  (s)G(s) 


-1 


Rearranging: 


-1 


G  (s)H(s)  =  G  (s)G(s)    -  I 

—  p    -»       —  p    -         - 

Then  multiplying  both  sides  of  equation  (16)  bv  G  (s) 

? 

H(s)  =  G(s)"1  -  G  (s)"1 


-1 


15) 


(17) 


(13) 


3 .   Cascade  Compensation 

Cascade  compensation  can  be  treated  in  a  similar  fashion 
Consider  the  same  basic  plant  with  cascade  compensation  as  shown 
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in  Figure   2-5. 


UKS] 


U2(S) 


CMS] 


C2(S) 


Figure  2-5 
Multivariable  Plant  with  Cascade  Comoensation 


The  algebraic  equations  for  R, (s)  and  R-(s)  are: 

Rl(s)  =  Gcll(s)U1(s)  +  Gci2(s)U2(s) 
R2(s)  =  Gc21(s)Ul(s)  +  Gc22(s)U2(s) 


or 


Rl(s) 
R2(s) 


Gcll(s) 
Gc21(s) 


Gcl2  (s 
Gc22(s 


Ui(s) 
U2(s) 


(19) 
(20) 

(21) 


in  the  general  case: 

C(s)   =   Gp(s)Gc(s)U(s) 
G(s)   =   Gp(s)Gc(s) 

If  Gp(s)  is  a  square  non-singular  matrix  the  compensation 
Gc(s)  is  given  by: 

Gc(s)  =  Gp(s) _1  G(s) 


(22) 
(23) 


(24) 
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The  compensated  plant  can  be  depicted  as  shown  in  Figure  2-6. 


R(S1 


U(S) 


CIS) 


> 


Figure  2-6 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Cascade  Compensation 


4.   Cascade  and  Feedback  Compensation. 

Cascade  and  feedback  compensation  can  also  be  combined 
as  shown  in  Figure  2-7. 


Figure  2-7 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Cascade  and  Feedback  Compensation 
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E(s)   -   R(s)  -  H(s)C(s)  (25) 

U(s)   =   Gc(s)E(s)  (26) 

C(s)   =   G  (s)U(s)   =   Gp(s)Gc(s)  [R(s)  -H(s)C(s)]  (27) 

C(s)   =   Gp(s)Gc(s)R(s)  -  Gp(s)Gc(s)H(s)C(s)         (28) 
Gp(s)Gc(s)R(s)   =   [I  +  Gp(s)Gc(s)H(s)]C(s)  (29) 

If  Gp(s)  is  an  ixj  matrix  Gc(s)  will  be  a  j  x  j  matrix,  and 

s 

H(s)  will  be  jxi,  in  which  case  Gp (s) Gc ( s) H (s)  is  a  square 
matrix  (ixi).   If  [I  +  "Gp  (s)  Gc  (s)  H-Cs)  ]  is  non-singular,  it  can  be 
inverted  and: 

C(s)   =   [I  +  Gp(s)Gc(s)H(s)  ]"1Gp(s)Gc(s)R  (s)        (30) 

C(s)   =   G(s)  R(s) 

The  compensated  system  transfer  matrix  G(s)  (ixj)  is: 

G(s)   =   [I  +  Gp(s)GC(s)H(s)  ]_1  Gp(s)Gc(s)  (31) 

Equation  (31)  contains  two  unknown  matrices,  H(s)  and  Gc(s) . 
One  of  these  must  be  selected  arbitrarily  and  then  in  certain 
special  cases  the  other  can  be  solved  explicitly.   For  example: 

h(s)  =  [i  !J  ■  j  (32) 

In  which  case: 


G(s)   =   [I  +  Gp(s)Gc(s) ]"X  Gp(s)Gc(s)  (33) 

[I  +  Gp(s)Gc(s)  ]G(s)   =   Gp(s)Gc(s)  (34) 

G(s)   =Gp(s)Gc(s)[I-G(s)]  (35) 

If  Gp(s)  is  a  square  matrix  G(s)  and  Gc(s)  will  be  square.   If 
Gp(s)  [G(s)~   -  I],  and  G(s)  are  non-singular: 
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Gp(s)"1   =  Gc(s)[G(s)"1  -  I]  (36) 

Gc(s)     -  Gp(s)"1  [G(s)"1  -  I]"1  (37) 

If  Gc(s)  and  {[G(s)    -  l]Gp(s)}  are  assumed  to  be  non-singular 
equation  (37)  can  be  rearranged  for  convenience: 

Gp(s)Gc(s)   =   [G(s)"1  -  I]"1  (38) 

Gp(s)   «   [G(s)"1  -  I]"1  Gc(s)"1  (39) 

Gc(s)"1   =   [G(s)"1  -  I]  Gp(s)  (40) 

Gc(s)   =   {[G(s)"1  -  I]  Gp(s)}"1  (41) 

This  result  can  also  be  obtained  directly  from  equation  (37) 
using  the  matrix  identity  (AB)~   =  B~   A~  . 

In  each  of  the  systems  discussed  above  the  compensation 
cannot  be  solved  for  in  closed  form,  except  in  special  cases 
as  indicated. 


B.   PARAMETER  OPTIMIZATION  -  COMPUTER  METHOD 

1.   Optimization  Technique 

The  Complex  Method  of  constrained  optimization  (BoxPLX) 
was  utilized  in  this  program.   The  algorithm   suggested  by 
Box  (9)  was  modified  by  the  programmer,  R.  R.  Hilleary,  Naval 
Postgraduate  School,  to  find  the  constrained  minimum,  and 
includes  an  integer  programming  option,  as  suggested  in  (10) . 
The  Complex  Method  is  a  modification  of  the  unconstrained,  direct 
search,  Simplex  Method  introduced  in  (11).   Direct  search  methods 
compare  function  and  constraint  values  only  in  searching  for 
the  minimum  value  of  the  function.   These  direct  search  methods 
have  been  widely  used  because  of  their  robustness,  reliability, 
and  ease  in  programming  and  use.   They  have  widespread  applic- 
ability.  The  one  drawback  is  that  they  are  generally  less 
efficient  than  gradiant-based  techniques  (12) .   A  more  detailed 
discussion  of  the  Complex  Method  is  included  in  the  program 
documentation. 

In  this  program,  the  main  program  simulates  the  reference 
output  and  then  calls  FUNCTION  BOXPLX  which  in  turn  calls  FUNCTION 
FE.   FUNCTION  FE  calls  SUBROUTINE  Plant  which  simulates  the 
compensated  system.   FE  then  computes  the  performance  measure 
and  returns  to  BOXPLX.   BOXPLX  keeps  track  of  the  compensator 
parameter  values,  computes  a  new  set  of  free  parameters  and 
calls  FE  again.   This  iteration  continues  until  termination  as 
described  in  the  documentation.   See  Appendix  A  for  more 
information  on  the  program. 
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2.   Performance  Measures 

In  order  for  a  numerical  method  for  parameter  opti- 
mization to  produce  a  result,  a  decision  process  must  be  com- 
pleted, the  result  of  which  is  the  set  of  parameters  which  is 
defined  as  "optimum."   This  term  optimum  is  an  enigima  of  sorts. 
The  program  addressed  here  utilizes  a  comparison  of  time  domain 
response  of  a  compensated  system  with  a  reference  response 
specified  by  the  user.   The  word  optimum  above  is  dependent  on 
the  designer's  knowledge  of  what  reference  response  is  best 
suited  to  the  output  he  is  trying  to  control.   Furthermore, 
with  the  exception  of  very  simple  linear  systems,  it  is  diffi- 
cult to  match  a  reference  response  exactly,  so  the  designer 
must  decide  if  the  comparison  at  steady  state  is  more  or  less 
important  than  that  during  the  transient  period.   Additionally 
very  high  frequencies  with  periods  less  than  one  half  the 
integration  step  size  may  be  filtered  out  by  the  digital  process 
in  the  computer,  and  not  be  present  in  the  simulated  output. 
The  designer  is  faced  with  some  decisions  which  he  must  make 
objectively  before  the  "optimum"  set  of  parameters  can  be 
computed. 

The  decision  process  mentioned  above  is  based  upon  a 
performance  measure  which  is  an  indicator  of  the  "goodness" 
of  a  given  set  of  parameters  being  varied  in  the  optimization 
process.   The  designer  must  select  the  performance  measure;  the 
performance  measure  is  then  minimized  by  numerical  methods , 
since  the  objective  here  is  to  have  the  system  response  match 
the  reference  response. 
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The  term  selectivity  is  used  to  describe  a  quality  of  the 
performance  measures.   If  three  dimensional  space  is  considered 
with  the  performance  measure  plotted  as  a  function  of  two  free 
variables  in  a  compensator  of  the  system,  the  resultant  topology 
should  include  one  or  more  minima  for  the  performance  measure. 
The  minima  with  the  lowest  value  are  referred  to  as  the  global 
minima,  although  the  minimization  performed  is  constrained. 
Selectivity  is  directly  related  to  the  steepness  of  the  slope 
of  the  contour  in  the  area  immediately  adjacent  to  the  minima. 
This  selectivity   is  a  function  of  the  performance  measure  used 
as  well  as  the  plant  and  compensator  dynamics  (2,3)  . 

Several  performance  measures  have  been  suggested  in  numerous 
sources,  some  of  which  will  be  discussed  here.   In  each  case 
the  performance  measure  will  be  represented  by  the  variable  J. 

Consider  a  single-input  single-output  system  with  output, 
c(t) ,  and  reference  response,  r(t) .   The  error,  e(t) ,  for  this 
case  is  defined  as  the  difference  between  the  reference  res- 
ponse and  the  system  response  or: 

e(t)   =   r(t)  -  c(t)  (42) 

Integral  square-error.   The  defining  equation  for  finite 
time  is : 

J   =        e2(t)dt  (43) 


■/: 


This  performance  measure  is  not  very  selective  and  tends  to 
emphasize  large  errors  and  ignore  small  ones  (2,8) . 
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Integral  of  time  multiplied  square-error.   The  defining 
equation  for  finite  time  is: 


-r 

Jo 


t  e2(t)dt  (44) 


Steady  state  error  is  penalized  more  heavily  by  this  per- 
formance measure.   This  performance  measure  is  said  to  be 
more  selective  than  the  integral  square-error  criterion  in  that 
the  value  of  the  integral  tends  to  change  more  rapidly  with 
changes  in  the  system  parameters  (2,8)  . 

Integral  absolute  error.   The  defining  equation  for  finite 
time  is : 


r, 

|e(t 

Jo 


J  =      |e(t) |dt  (45) 

Jo 

This  performance  measure  is  slightly  more  selective  than  the 

integral  square-error  criterion  (2,3). 

Integral  of  time-multiplied  absolute  error.   The  defining 

equation  for  finite  time  is : 

r 

J   =  I     t] e(t) |dt  (46) 

This  performance  measure  is  more  selective  than  those  mentioned 
previously  and,  as  in  the  case  with  the  integral  of  time  multi- 
plied square-error  criterion,  this  performance  measure  tends  to 
emphasize  errors  late  in  the  transient  period  and  in  steady 
state  and  deemphasize  those  which  occur  early  in  the  transient 
period  (2,3) . 
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The  optimization  of  the  performance  measure  can  be  addressed 
in  two  categories;  the  case  where  the  reference  response  speci- 
fied is  the  forcing  function;  or  the  reference  response  is  some 
function,  usually  a  second  order  response,  which  the  designer 
believes  to  be  a  model  response  for  his  application.   In  the 
first  case  the  peak  overshoot,  setting  time,  rise  time,  and 
natural  frequency,  are  determined  by  the  dynamics  of  the  com- 
pensated system  and  by  the  performance  measure  selected.   How- 
ever, if  the  designer  desires  a  second  order  response  with  a 
certain  minimal  rise  time,  he  can  specify  a  second  order 
reference  response  which  meets  this  requirement,  and  use  the 
integral  absolute  error  or  integral  error-squared  performance 
criteria  to  optimize  the  compensation  parameters . 

Provisions  for  both  of  these  approaches  have  been  included 
in  the  program  in  that  the  user  can  reference  a  second  order 
response,  where  he  specifies  the  damping  factor,  natural  frequency 
and  gain,  or  he  can  utilize  the  table  loop-up  feature  in  which 
he  can  specify  any  reference  response  he  desires. 

Other  performance  measures  can  be  developed  using  classical 
measures  of  system  performance  such  as  maximum  overshoot,  time 
for  the  error  to  reach  its  first  zero,  time  to  reach  maximum 
overshoot,  settling  time,  frequency  of  oscillation  of  the 
transient  or  steady  state  error.   However,  performance  measures 
based  on  these  characteristics  would  be  complex  and  increase 
computation  time,  and  their  desirability  over  those  previously 
discussed  is  questionable. 

The  performance  measure  used  in  t*he  optimization  program 


(see  Appendix)  is  the  sum  of  the  weighted  performance  measures 
of  each  output  of  the  multivariable  system. 

M 
J   =    Z  (Wi) (Jik)  (47) 

i=l 

The  subscript  i  denotes  the  output  for  which  J.  .   is  calculated 

1  z  K 

The  subscript  k  denotes  one  of  four  performance  measures  which 

can  be  selected  by  the  user  to  be  applied  to  each  output  in 

turn : 

Integral  square  error 

N 

Ji,l   =    S    (Rj,i  -  C.  .) 2  (48) 

j-i  D'L 

Integral  of  time  multiplied  square  error 

N 

Ji,2   =    Z         ((Rj,i  -  C.  .)2)(iAT)  (49) 

j=l  3'1 


Integral  of  absolute-error 

N 

Ji,3   -   Z      |  Rj,i  -  C.  .  |  (50) 

j=l  ^'L 

Integral  of  time  multiplied  square-error 

N 

Ji/4   =    Z       |  Rj,i  -  C.  . |  (iAT)  (51) 

J  t  ^ 

j-l 

N  represents  the  number  of  integration  steps  in  the  simulation, 
AT  represents  the  integration  step  size,  M  is  the  number  of 
outputs  being  optimized,  and  W  is  the  weighting  factor  applied 
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to  each  output. 

It  should  be  noted  that  the  performance  measures  are  an 
approximation  and  when  the  weighting  factor  is  introduced 
the  performance  measure  diverges  significantly  from  the  defining 
integral  equation.   Comparison  of  performance  measures  for  a 
given  plant  with  various  types  of  compensation  should  be  done 
with  the  same  weighting  factors.   The  same  performance  measure 
could  be  achieved  for  a  system  using  two  different  compensation 
schemes,  since  the  performance  measures  for  each  channel  may 
vary  significantly. 
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C.   A  TWO  INPUT-TWO  OUTPUT  SYSTEM  EXAMPLE 

1.   Transfer  Matrix  Solution 

In  order  to  demonstrate  the  program  it  was  decided  to 
analyze  a  simple  system,  synthesize  the  compensation  required 
to  produce  a  specified  output,  and  then  show  that  the  optimi- 
zation program  will  arrive  at  the  same  results.   It  is  emphasized 
that  this  program  will  not  determine  the  type  of  compensation 
required  to  achieve  the  desired  response.   It  will  set  free 
parameters  in  the  compensators  (specified  by  the  user)  to  most 
nearly  produce  the  desired  response. 

The  analysis  and  compensation  was  performed  using  the 
transfer  matrix  method  described  by  OGATA  (2) .   The  uncompensated 
plant  and  the  compensated  plant  were  simulated  to  show  that  the 
compensation  achieved  the  design  specifications.   Certain  para- 
meters in  the  compensator  were  defined  as  free  parameters. 
These  free  parameters  were  offset  from  their  optimal  values 
and  constraints  were  introduced  which  restricted  the  range  of 
values  for  these  parameters  during  optimization.   The  plant 
and  compensators  with  free  parameters  were  simulated  in  the 
optimization  program  and  optimization  was  accomplished. 

To  illustrate,  a  simple  two-input  two-output  linear 
system  shown  in  Figure  3-1  was  selected  as  the  plant  to  be 
compensated.   The  design  specifications  chosen  were: 

1.  Channel  one  and  channel  two  are  to  be  decoupled  during 
the  transient  period  as  well  as  in  steady  state. 

2.  Channel  one  is  to  have  a  second  order  response  to 

a  step  input  with  a  natural  frequency,  W  ,  of  10,  a 
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Figure  3-1   Multivariable  Plant  to  be  Compensated 


Figure  3-2   Multivariable  Plant  with  Generalized  Compensation 
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damping  factor,  5 ,  of  0.4,  and  no  steady  state  error. 

3.   Channel  two  is  to  have  a  second  order  response  to  a 
step  input  with  a  natural  frequency,  W  ,  of  4,  a  damping 
factor,  z, ,  of  0.6,  and  no  steady  state  error. 

To  achieve  total  decoupling  the  off-diagonal  elements  of 
the  transfer  matrix  of  the  compensated  system  must  be  zero. 
The  desired  compensated  system  transfer  matrix,  G(s),  is: 


G(s) 


100 


S   +  8S  +  100 


0 


16 


5   +  4  .  8S 


16 


(52) 


The  input/output  equations  obtained  from  Figure  3-1  are: 


y1(s: 


Ul(s)   (S+7) (S+12) 


+  Vs)   UTTIT 


(53) 


Y2(3)   =   U2(S)   (s+2) (s+4.3 


Yx(s 


(s  +  2) 


(54) 


Substitution  of  equation  (53)  into  equation  (54),  and  equation 
(54)  into  equation  (53)  produces  equations  (55)  and  (56) . 


Y1(s) 


U]_(s) 


+ 


2U2(s) 


+ 


3Y1(s) 


s+7)  (s+12)    (s  +  2)  (s+4.3)  (s  +  12)    (s  +  2)  (s  +  12) 


(55) 


a,  (a) 

y  (s)  =  : 

*2^;    (s+2)  (s  +  4.3) 


4Ux(s 


8Y2  (s) 


s+12) (s+7) (s+2)    (s+2) (s+12) 


(56) 


Collecting  like  terms  in  equations  (55)  and  (56)  and  rearranging 
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s^+14s+16 
lls;  (s+2) (s+12) 


U1(s) 


2U2(s) 


(s+7)  (s+12)    (s  +  2)  (s  +  4.3)  (s+12) 


(57) 


4Ux(s) 


U2(s) 


t    \    s  +14s+16     

*2{S)     (s+2) (s+12)  '   (s+12) (s+7) (s+2)  T  (s+2) (s+4.3) 


(58) 


Multiplying  both  sides  of  equations  (57)  and  (58)  by 
(s+2) (s+12) 


s  +14s+16 


Yx(s)  = 


(s+2)U1(s) 


(s+7) (s  +14s+16 


+ 


2U2(s) 


(s+4.3) (s  +14s+16) 


(59) 


Y2(s)  = 


4U. (s) 
I 


(s  +  12)U2(s) 


(s+7) (s2+14s+16)    (s+4.3) (s2+14s+16) 


(60) 


The  transfer  Matrix  for  the  uncompensated  plant  is: 


s+2 


G  (.) 


(s+7) (s  +14s+16)    (s+4.3) (s^+14s+16) 


(s+12 


(s+7) (s  +14s+16)    (s+4.3) (2^+14s+16 
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The  generalized  compensation  to  be  used  is  shown  in  Figure 
3-2.   Equation  (41)  will  be  used  to  calculate  the  required 
compensation;  since  in  this  case  the  feedback  compensation 
transfer  matrix,  H(s),  is  the  identity  matrix. 


Gc(s)   =   {[G(s)""1  -  ?]9p(s)  >~X 


(41) 


Since  G  (s^  is  a  2x2  matrix  in  this  example  G  (s)  and  G(s)  will 
~p  r         -.c 

also  be  2x2.   By  inspection  (Equation  (52)),  G(s)  is  non-singular 

and  can  be  inverted. 
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The  compensated  system  transfer  matrix  G(s)  is  given  in 
equation  (52).   The  inverse  of  G(s)  is  calculated  by  equation 
(62)  . 


G(s) 


-1       COF' 


G(s) 


(62) 


Where  COF   is  the  transpose  of  the  cof actor  matrix  of  G(s), 
and  |G(s)|  is  the  determinant  of  G(s). 

16  0 


COF' 


s  +4.8s+16 


100 


s  +8s+100 


(63) 


G(s)  |  = 


(16) (100) 


(s2+8s+100) (s2+4.8s+16) 


(64) 


G  is) 


-1 


s  ^-8s  +  100 
10  0 


5  +4.35+16 

16 


(65) 


[G(s)"1-!]   = 


s  +8s 
100 


L 


s  +4 . 8s 
16 


(66) 


The  plant  transfer  matrix,  G  (s),  is  given  in  equation  (61)  as 
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V~s)  - 


(s+2) 2 

(s+7) (s2+14s+16)    (s+4.3) (s2+14s+16) 


(s+12) 


(s+7) (s2+14s+16)    (s+4.3) (s2+14s+16 


(61) 


Multiplying  equation  (61)  by  equation  (66) : 


[  GtsO^-IjG  (s) 


s  (s+8)  (s  +  2) 


2s(s  +  8) 


100(s+7) (s2+14s+16)   100(s+4.3) (s2+14s+16) 


(67) 


4(s+4.8) 


s(s+4.8) (s+12) 


16 (s+7) (s2+14s+16)   16(s+4.3) (s2+14s+16) 


If  this  matrix  is  non-singular  it  can  be  inverted 


_i 


[G(s)   -I]Gn(s) 


s  (s+3)  (S+2)  (s+4.3)  (s  +  12) 

(16) (100) (s+7) (s+4.3) (s2+14s+16)2 


(63) 


3a     (s+4.8) (s+8 


(16)  (100)  (s+7)  (s  +  4.3)  (s2  +  I4sfl5)2 


[g(s: 


-i 


I]G    (a) 

-    -P    - 


s2(s+4)  (s  +  3)  (s2+14s+24-3) 

(16) (100) (s+7) (s+4.3) (s2+14s+16)2 


(69) 


s     (s+4.8) (s+3) 


(16)  (100)  (s+7)  (s+4.3)  (s2  +  14s  +  16) 
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r-    s(s+4.8)  (s+12) 

~~2 


COF 


2s(s+8) 


16(s+4-3)(s+14s+16)    100(s+4.3)(s+14s+16) 


-4s(s+4.8) 


s  (s+8)  (s+2) 


_  16(s+7) (s  +14s+16) 


100(s+7) (s  +14s+16)   J 


(70) 


9c  ^ 


{     [G(s)"1  -  l]G  (a)}"1 


100 (s+7) (s+12) 
s (s+8) 


-400(s+4.3) 


s(s+8) 


-32(s+7) 
s(s+4.8) 


16(s+2) (s+4.3) 
s(s+4.3) 


(71) 


One  of  the  possible  realizations  of  this  compensation  is 
shown  in  Figure  3-3.   To  obtain  a  confirmation  of  the  design, 
the  system  was  simulated  with  and  without  compensation  using 
International  Business  Machines'  Digital  Simulation  Language, 
DSL  (13).   Figures  3-4  through  3-9  show  the  results  of  that 
simulation,  which  indicate  that  the  compensators  used  do  in 
fact  achieve  the  stated  objectives  of  the  design.   In  each  plot, 
trace  number  1  is  the  desired  response,  Z,  and  trace  number  2 
the  system  response,  Y.   In  Figures  3-4B,  3-73,  3-83,  and  3-93 
the  two  traces  are  superimposed.   In  Figures  3-5A  and  3-6A 
the  desired  responses  are  zero.   The  reader  should  note  the 
scale  factors  present  on  some  of  these  plots. 

2.   Parameter  Optimization  Solution 

The  compensated  system  and  the  desired  response  curves 
were  then  simulated  and  solutions  were  obtained  using  the 
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Figure  3-4A 
Yl  and  Zl  vs  Time  with  Rl  =  Unit  Step,  R2 
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Figure    3-4B 

Yl  and  Zl  vs  Time  with  Rl  =  Unit  Step,  R2 
For  Compensated  System 
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Figure  3-5A 

Y2  and  Z2  vs  Time  with  Rl  =  Unit  Step,  R2  =  0.0 

(Y2    =  3.57  X  10"2) 
ss 


~— 


o.oo 


2.00 


4.00 


S.OO 


8.00 


10.00 


Figure  3-53 

Y2  and  Z2  vs  Time  with  Rl  =  Unit  Step,  R2 
For  Compensated  System 
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Figure    3-6A 
Yl   and    Zl   vs   Time   with    Rl   =    0.0,    R2    =   Unit   Step 

(Yl        =    2.91    X    10"2) 
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Figure    3-63 

Yl  and  Zl  vs  Time  with  Rl  =  0.0,  R2  =  Unit  Step 

For  Compensated  System 
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Figure  3-7A 
Y2  and  £2  vs  Time  with  Rl  =  0.0,  R2  =  Unit  Step 
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Figure    3-73 

Y2  and  Z2  vs  Time  with  Rl  =  0.0,  R2  =  Unit  Step 

For  Compensated  System 
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Figure    3-3A 
Yl   and    Zl   vs   Time  with   Rl   =   Unit   Step,    R2    =   Unit   Step 
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Figure    3-3B 

Yl  and  Zl  vs  Time  with  Rl  =  Unit  Step,  R2  =  Unit  Step 

For  Compensated  System 
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Figure    3-9A 
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Y2  and  Z2  vs  Time  with  Rl  =  Unit  Step,  R2  =  Unit  Step 
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Figure  3-9B 

Y2  and  Z2  vs  Time  with  Rl  =  Unit  Step,  R2  =  Unit  Step 

For  Compensated  System 
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program  discussed  in  this  section.   The  time  response  results 
are  shown  in  Figure  3-10.   There  is  a  slight  difference  between 
the  desired  response  curve  and  the  output  of  the  simulated 
system. 

Next,  two  parameter  optimization  was  performed.   The  pro- 
portional gain  constant  of  each  compensator  was  selected  as 
the  adjustable  parameters.   Each  of  the  four  performance  measures 
given  by  equations  (47)  through  (51)  were  used  in  successive 
optimizations.   All  other  aspects  of  the  program  were  unchanged 
during  this  optimization.   Equal  weighting  was  applied  to  each 
channel.   The  gains  10  0  and  16  were  chosen  because  their 
variation  should  have  similar  impact  on  the  two  respective  out- 
puts.  The  upper  limits  were  10  5  and  21  and  the  lower  limits 
were  95  and  11  respectively.   The  integration  step  size  was 
.02  and  the  final  time  was  5  seconds.   Optimization  was  begun 
with  the  two  gains  set  at  10  5  and  11  respectively,  and  Rl  =  R2  = 
Unit  Step  functions.   The  results  are  given  in  Table  3-1. 

The  system  was  simulated  using  these  values  for  the  two 
free  parameters.   The  simulation  results  are  shown  in  Figures 
3-11  through  3-14.   This  simulation  indicates  the  time  multiplied 
integral  square-error  performance  measure  produced  the  best 
results . 
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Value  of 
Performance        Minimum  G100  G16 

Measure  Performance 

Measure 


Integral  -> 

Square-  1.2026097  x  10       104.9999      20.99998 

Error 


Time 

Multiplied  - 

Integral  1.641641   x  10       102.2554      15.27091 

Square- 
Error 


Integral  , 

Absolute  1.353197   x  10       104.9862      14.32275 

Error 


Time 

Multiplied  . 

Integral  9.994698   x  10       103.2723      14.31198 

Absolute 

Error 


Table  3-1 


RESULTS  OF  OPTIMIZATION 
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Figure    3-10A 


Yl  and  Reference  1  vs  Time  Simulated  with  Thesis  Program 
Parameters  Set  at  Analytical  Optimum 
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Figure    3-10B 

Y2  and  Reference  2  vs  Time  Simulated  with  Thesis  Program 
Parameters  Set  at  Analytical  Optimum 
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Figure  3-11A 

Yl  and  Reference  1  vs  Time  After  Optimization  With 
Integral-Square  Error  Performance  Measure 
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Figure  3-113 

Y2  and  Reference  2  vs  Time  After  Optimization  With 
Integral  Square-Error  Performance  Measure 
(Curve  With  Greatest  Overshoot  is  ZZ) 
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Figure  3-12A 

Yl  and  Reference  1  vs  Time  After  Optimization  With 
Time  Multiplied  Integral-Square  Error  Performance  Measure 
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Figure  3-12B 

Y2  and  Reference  2  vs  Time  After  Optimization  With 
Time  Multiplied  Integral  Square-Error  Performance  Measure 
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Figure  3-13A 

Yl  and  Reference  1  vs  Time  After  Optimization  With 
Integral  Absolute  Error  Performance  Measure 
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Figure    3-13B 

Y2  and  Reference  2  vs  Time  After  Optimization  With 
Integral  Absolute  Error  Performance  Measure 
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Figure  3-14A 
Yl  and  Reference  1  vs  Time  After  Optimization  With 


Time  Multiplied  Integral  Absolute  Error  Performance 


Measure 
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Figure  3-143 
Y2  and  Reference  2  vs  Time  After  Optimization  With 


Time  Multiplied  Integral  Absolute  Error  Performance 


Measure 
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CONCLUSIONS 

From  the  results  of  the  two-input  two-output  example 
system,  the  following  conclusions  can  be  stated: 

1.  The  technique  of  cost  function  minimization  in  the 
time  domain  can  be  used  effectively  for  compensator  parameter 
optimization  in  multiple  input  multiple  output  control  systems 

2.  Cost  function  minimization  is  not  a  substitute  for 
control  system  design.   The  user  must  be  able  to  select  the 
proper  type  of  compensator  which  is  capable  of  achieving  the 
desired  results,  before  initiating  the  optimization  of  the 
free  parameters . 

3.  The  four  performance  measures  evaluated  produced 
slightly  different  solutions  to  the  same  problem,  with  the 
time  multiplied  integral  square-error  yielding  the  best 
results  in  the  example  given. 
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APPENDIX   A 

£*********************************«**********  ***  ****  *****  v*C 

c  c 

C  COMPENSATOR    OPTIMIZATION    IN                                             C 

C  C 

C  MULTIVARIA3LE    CONTROL    SYSTEMS                                          C 

c  c 

C  3Y                                                                             C 

c  c 

C  JOHN    T.    MOWREY                                                             C 

c  c 

C  MAJOR    USMC                                                                C 

C  C 

C  DECEMBER    1978                                                             C 

C  C 

C  THIS  PROGRAM  WAS  DEVELOPED  AS  A  GENERAL  PURPOSE  DESIGN  C 
C  AID  TO  9E  USED  TO  OPTIMIZE  FREE  PARAMETERS  IN  THE  COMP-  C 
C    ENSATCRS    OF    A    MLLTI VARI ABLE    CONTROL    SYSTEM.     IT    CAN    ALSO       C 

C  BE  USED  FTP  THE  OPTIMIZATION  OF' FREE  PARAMETERS  IN  A  C 
C    SINGLE    INPUT    SINGLE   OUTPUT    SYSTEM'  OR    xO    SIMULATE    EITHER       C 

C  TYPE  OF  SYSTEM.  OPTIMIZATION  IS  ACCOMPLISHED  BY  MINIMIZ-  C 
C    IMG    THE    ERROR     BETWEEN    A    USER    SPECIFIED    REFERENCE    CURVE         C 

C    A  NO    THE    OUTPUT    OF    THE    SIMULATED    COMPENSATED    SYSTEM.  C 

C  C 

Q  *****************************************************  *****(; 

c  c 

C  DATA    CARDS                                                                C 

C  **2C*********                                                                                                                                         £ 

C  C 

C    CARD  1      MRUNStNOPT,NGRAPrifNREF,IINTNOLTfIPRINTfIFREC               C 

C  C 

C  FORMAT(SI5)                                                                                                         C 

c  c 

C  NRUNS-FLAG    FOR    NUMBER    OF    RUNS.    NRUNS    MUST    EQUAL       C 

C  1     WHEN    OPTIMI  ZING.                                                                                         C 

C  C 

C  NOPT-FLAG    FOR    OPTIMIZATION/SIMULATION                                C 

C  1-OPTIMIZAriON                                                                                       C 

C  2-SIMULATION    ONLY                                                                                  C 

C  (IF    NOPT  =  l    THE     PLANT    WILL     BE    SIMULATED    AFTER         C 

C  OPTIMIZATION    AND    THE    REQUESTED    OUTPUT    PROCUCEC    C 

c  c 

C  NGRAPH-FLAG    FOR    TY°E    OF    GRAPHICAL    OUTPUT    DESIRED    C 

C  1-PRODUCES    PRINTER    =>LQTS                                                              C 

C  2-oc>0  DUCES    VERSA  TEC     PL3TS                                                             C 

C  3 -MO    GRAPHS     PRODUCED                                                                          C 

c  c 

C  NPEF-FLAG    FOR    REFERENCE    CURVE                                                     C 

C  1 -PRO  DUCES     CURVE                                                                                  C 

C  2-NO    REFERENCE    CURVE     IS    PRODUCED                                          C 

C  C 

C  IIN-THE    NUMBER    OF    INPUTS     (UP    tq    25)                                      C 

^  r 

C  NOUT-THE    NUMBER    OF    OUTPUTS     (U°    TO    3)                                     C 

c  c 

C  IPRINT-FLAG    FOR    TABULATED    DATA    FOR    REFERENCE               C 

C  CURVE    AND   SYSTEM    RESPONSE.                                                                C 

C  1-PRODUCES    PRINTED    OUTPUT                                                             C 

C  2-NO    OUTPUT    DA^A    PRINTED                                                                C 

c  c 

C  I5=REQ-CREQUENCY    OF    PRINTED    OUTPUT.     IF     [FREO-10          C 

C  EVE^Y    10TH    DATA     POINT    ^1 LL     BE    PRINTED.                          C 

C  C 

c  c 

C    CARD  2      T47,DT,TF                                                                                                              C 

c  c 

C  FORMAT(3F10»5>                                                                                                 C 

C  C 

C  T4-7-INITIAL    TIME                                                                                          C 

c  c 

C  DT-INTEGRATION    STEP    SIZE 

C  r 

C  TF-FINAL    TIME                                                                                                  C 
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c 


c  c 

c  c 

C    CARO    3       NV,NAV  tNTAfNPRt  IPtlPM  C 

c  c 

C  F3RMAT(6I5)                                                                                                         C 

C  C 

C  THIS   CARD    MUST   3E    CHITTED    IF    N0PT=2                                      C 

c  c 

C  NV-THE    NUMBER    OF    FREE   VARIABLES    T3    BE   SET    BY              C 

C  OPTIMIZATION.                                                                                           C 

C  C 

C  NAV-THE    NUMBER    OF    AUXILLIARY    VARIABLES                               C 

C  C 

C  NTA-THE    NUMBER   OF    TRIALS    ALLOWED    30XPLX                           C 

C  C 

C  NPR-FREQUENCY    OF    OUTPUT    FROM    BOXPLX    FOR                            C 

C  DIAGNOSTIC    PURPOSES                                                                          C 

c  c 

C  IP-FLAG    FOR    INTEGER /REAL*8    OPTIMIZATION.                         C 

C  C 

c  *               *                                           c 

C  *       NOTE       *                                                                C 

C  *                      *                                                                C 

r 

SEE  BOXPLX  FOR  ADDITIONAL  INFORMATION  ON  THE  5  VARIABLES  C 
C  LISTEO  ABOVE.  30XPLX  WAS  CONVERTED  TO  DOUBLE  PRECISION  C 
C    ^OR    USE    IN     THIS    PROGRAM;     HOWEVER    THE    DOCUMENTATION    WAS  C 

C    NOT    CHANGED  C 

C  C 

C  IPM-R.AG    FOR    PERFORMANCE    MEASURE.                                           C 

C  C 

C  THE    PERFORMANCE    MEASURES    USED    ARE    PROPORTIONAL    C 

C  TO    THOSE    INDICATED    BELOW.                                                             C 

C  C 

C  1-INTEGRAL    SQUARE-ERROR     (WEIGHTED)                                   C 

C  2-TIME    MULTIPLIED    INTEGRAL    SCUARE-ERROR                      C 

C  (WEIGHTED)                                                                                             C 

C  3-INTEGRAL    ABSCLUTE    ERROR     (WEIGHTED)                              C 

C  4-TIME    MULTIPLIED     INTEGRAL    ABSOLUTE    ERROR                 C 

C  (WEIGHTED)                                                                                             C 

C  '            C 

C  C 

C    CARDS    4—7    XS(I  )  C 

c  ■** 

C  FORM-AT(8Fi0.5)                                                                                                 C 

c  c 

C  ALL    OF    THESE    CARDS    mijst    5=    OMITTED    Ic    NOPT=2              C 

C  XSdl     ARE    THE     STARTING    VALUES    3F    THE    FREE    PARA-      C 

C  METERS.    THE    NUMBER    OF    CARDS    REQUIRED    IS    OEPEN-          C 

C  DENT    ON    NV.     ONE    VALUE    MUST     BE    REAO    IN    FCR    EACH          C 

C  FREE    PARAMETER.     IE     IF    NV= 8    ONLY    CARO    4     IS    RE-            C 

C  QUIRED    AND    3     PARAMETER    VALUES    SHOULD    3E    ON     IT            C 

C  IN     8F1Q.  5    FDR  MAT.                                                                                            C 

C  C 

C  C 

C    CARDS     8-11     XU(  I  )  C 

C  C 

C  FORMAT     (8F10.5)                                                                                              C 

C  C 

C  THESE    CARDS    }RE    THE     SAME    AS    CARDS    4-7    EXCEPT               C 

C  XU(I)     ARE   THE    iJPOER    LIMITS    ON    THE    FREE    PARA-               C 

C  METERS.    OMIT    IF    NQPT*2.                                                                      C 

C  C 

C  C 

C    CARD  S    12-15    XL(  I)  C 

C  C 

C  FORM/ST(3FI0.5  )                                                                                                C 

C  C 

C  THESE    CARDS    ARE    THE    SAME    AS    CARDS    4-7    EXCE^               C 

C  XL(I)    ARE    THE    LOWER    LIMITS    ON    THE    FREE     PA  R  A- 

C  METERS.     OMIT    IF     NCPT    =2.                                                                      C 

C  C 
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C    CAROS    16,18,     AND    20     131)7,  I WT,  I  TAB  C 

C  C 

C  =  QRMATt 3F10.5)  C 

C  C 

C  I  OUT-THE    NUMBER    OF    THE    BLOCK     FROM    WHICH   THE    OUT-  C 

C  °UT    IS    TAKEN,     I F    2    OR    MORE    BLOCKS    FEED    AN    OUT-  C 

C  PUT    NODE, A    TY°E    1    BLOCK     WITH    G=l    MUST    BE  C 

C  PLACED    BETWEEN    THAT    NODE    AND    THE    OUTPUT.  C 

C  C 

C  FOR    A   3    OUTPUT    SYSTEM    ALL    3    CARDS    ARE    REQUIRED.  C 

C  FOR    A    2    OUTPUT    SYSTEM    IS     AND    13    ARE    REQUIRED.  C 

r  Q 

C  IWT-AN    INTEGER    WEIGHTING    FUNCTION    WHICH    WILL  C 

C  BE    CONVERTED    TO    RgAL*8.AND    WHICH    ALLOWS    THE  C 

C  USE*    TO    PENALIZE    THE    ERROR    BETWEEN    THE    SYSTEM  C 

C  OUTPUT    AND    THE    REFERENCE    RES°0NSE    MORE    HEAVILY  C 

C  AT    ONE    OUTPUT    THAN   ANOTHER.    INTEGERS    BETWEEN     1  C 

C  AND    10    ARE    RECOMMENDED.     WEIGHTING    SHOULD    BE  C 

C  REDUCED    IF    OVERFLOW    IS     ENCOUNTERED.  C 

C  C 

C  ITAB-A    FLAG    FOR    TABLE    LOOK    UP    OF    THE    REFERENCE  C 

C  CURVFS.  C 

C  1-DROGRAM    WILL    READ    TABULATED    DATA    FOR    'HE  C 

C  REFERENCE    CURVES  .  C 

C  2-PROGRAM    WILL    CALCULATE    A    SECGNC    ORDER    RES-  C 

C  =>ONSE    CURVE.  C 

C  c 

c  c 

C    CARDS    17,19,     AND    21-3  ET  A  ,  DELTA  ,  WN  ,  AM°  ,QEL  AY  ,  I  NP  UT  C 

C  C 

C  C0RMAT{5F1Q.5,A4.)  C 

c  c 

C  THESE    CARDS     INPUT    THE    DA^A     REGUIRED     TO    CQMO|jTE    A  C 

C  SECOND    QRDE*    RESPONSE    C'JQVE.     1=    THE    USER    DESIRES  C 

C  THE    TABLE    LOOK   UP    FEATURE     THESE    CARDS    ARE    RE-  C 

C  PLACED    BY    TA8ULATED    DATA.     IF    A    TABULATED    REF-  C 

C  ERENCE    CURVE    IS    DESIRED    FOR    OLTPUT    NC.    1,     CARD  C 

C  17    IS    REPLACED    BY    TABULATED    DATA    IN    cF 1 0. 5    FORMATC 

C  THE    NUMBER    3F    DATA    POINTS     REQUIRED    IS    T  F/  DT     (TF  C 

C  AND    DT    FROM     CARD    2) .  C 

C  C 

C  C(I)                                   BcTA(I)  C 

C  R(I)            (S**Z+2*0ELTA< I )*WN{ I)*S+WN(I )**2)  C 

C  .  C 

C  CU)     IS    THE    DESIRED    REFERENCE    RESPONSE    ASSOC-  C 

C  IATED    WITH    OUTPUT    IOUT(I).  C 

C  P.U)     IS    THE    FORCING    FUNCTION    S  =>EC  I F I  EC    BY  :  C 

C  AM°(I)-THE    AMPLITUDE    OF     THE    FORCING    FUNCTION  C 

C  DELAYd  )-DELAY    A  F^  ER    ~M    BEFORE     THE    FORCING  C 

C  FUNCTION    IS     APPLIED(TiME    DOMAIN).  C 

C  INPUT(I)-THE    TYPE    OF    FORCING    F  UNCTION  :  S  TE=>  ,  C 

C  RAMP,     OR    3ARA.  C 

C  3  ETA (I) -THE    GAIN    OF    THE    TPANSFER    FUNCTION.  C 

C  OELTA(I)-THE    DAMPING    FACTOR.  C 

C  WN(I  )-THE   NATURAL    FREQUENCY.  C 

C  C 

C    CARD    22    M  C 

C  C 

C  F0RMAT(I2J  C 

C  C 

C  THE    NUMBER    OF    BLOCKS     IN    THE    SYSTEM     (UP    TO     23).  C 

C    CARDS    23-47    IC(  I )  f  ID(  I  ),  Ic(  I  ),  IV(  I  )  ,G  (I)  t  *  i  I )  ,  Z<  I)  C 


lw 


r 


:  c 

C  FORMAT(  412,  2X,3F1C5)                                                                              C 

C  c 

C  THE     NUMBER    3F     CAROS    REQUIRED     IS    DETE3*INED    BY    N.     C 
C 


^ 


C  IC(I)-THE    NUMBER    OF    THE    3L0CK(l-25).  C 

Z  C 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

r* 
w 

c 
c 
c 

r 

c 

r 

w 

c 
c 

c 
c 
c 

r 

c 

c 

r 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

r 

c 

c 
c 
c 
c 
c 
c 


IO(I)-THE    TYPE   OF    BLOCK,    S  I X    TYPES    ARE    AVAILABLE 

ic(i)     iv  m 

I  D(  II  *!  I         G  I 


ID(I)=2 


10(1)  =3 


IE(  I) 


IE<  I) 


G  I     IV(  I  ) 

S+P  I 


S*Z        I     IV(  I  ) 

G 

s+° 


IE(I  ) 
I  0  ( I  )  =4 


5**2+2*DELTA*fcN*S+WN**2 


I  V(  I  ) 


10(1 ) -5 


IE(I)       I 


U  L       / 

Mi 

/ 
/ 

/       LL 


I  V(  I) 


IV(I)=G*IE(I)    FOR  G*IE<I)  LE    LL 

ANO  G*IE( I )  GE    LL 

IV(  I)  =UL  FOR  G*IEU  )  GT    UL 

IV(  I)  =  LL  FOR  G*IE(  I)  LT    UL 


10(1) =6 


IE(  I) 


/  I 
/     I 

B    /        I 
I 


IV(  I) 


/    -8 
/  jG 
/- 


IV(I)aG*IE(I)    C0R    G*IE(I)     GE    6 

OR    LE    -3 
IV(I)=0     FOR    G*IE(I)    LT    3     ANO    GT    -3 

IE(I)-THE     INPUT    NODE    NUMBE3. 

IV(I)-THE    DUTPUT    NODE    NUMBER. 

G(I)-THE    GAIN     FIELu,    3(1)    REPRESENTS    THE    GAIN 
OF    THE    3L0CK. 

P(  I  )-THE    POLE    FIELD, 

FOR    I0(  I  )  =2    OP    3     P(I)sPGLE 
FnR    ID(I)=<V    P(  I  )*OAMPI.NG    FACTOR 
FOB    10(1)35     P<n  =  UP*ER    LI*IT,UL. 
=  0R    I0(I)=6    P(  I  PREFERENCE,  8. 

Z( I l-THE    ZERO  FIELO 

FOR    10(1 ) =3  Z(  I  ) =ZERO 

FOR    10(1)=-+  MIMNATUPAL    FP  EQUENCY  »WN. 

FOR    10(1) =5  Z(I)=LCWER    LIMIT,LL. 


C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
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c 
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c 
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-C 
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r 
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c  c 

C    CAROS    48-72    AMP ( I ) , DEL AY(  I ) , INPUT(  I ) f IEE ( I )  C 

c  c 

C  FORMAT(2F10.5 fA4,  12  )  C 

C  C 

C  ONE    CARO    IS    NEEDED    FOR    EACH    FCRCING    FUNCTION    ON  C 

C  THE    SIMULATED    COMPENSATED    SYSTEM.  C 

C  C 

C  AMP(  I), DELAY( I), ANDINPUT(  I  )     APE    AS    DESCRIBED    FOR  C 

C  CARDS    17,i9,ANC    21, EXCEPT    THEV    APPLY    TC    THE    COM-  C 

C  PENSATED    SYSTEM    IN    THIS    CASE.  C 

C  C 

C  ISE(I)-THE    NODE    TC    WHICH    THE    FORCING    FUNCTION     IS  C 

C  APPLIED.  C 

C  C 

C    CARDS    73-78    TITLE   CARDS    FOR    VERSATEC    PLCTS.     IF    NGRAPH<,  C 

C  NOUT*2    TITLE    CARDS    ARE    REQUIRED(2    PER    OUTPUT)  C 

C  TITLES    GO    IN    COLUMNS   1-4-8    ONLY.  C 

C  '  c 

c  *  *  c 

C  *       NOTE       *  C 

c  *  *  c 

r  r 

C  WHEN  CPTI  mi  zing  tONE  FORTRAN  CARD  °ER  FREF  VARIABLE  MUST  C 
C  BE  INSERTED  IN  SUBROUTINE  PLANT.  THE  APPROPRIATE  POSITIONC 
C  IS  IDENTIFIED  BY  THE  COMMENT  CARD  •  ENTEP  G< I )*C<  I  I  VAL-  C 
C    U5S   AT   THIS    POINT    ••  C 

c  c 

c  c 

r  M„»aji±3ji  ±*******  **************  *z*****i***^  ********  ****$**£ 
C  C 

C    "HE    FOLLOWING    J  CL    CARCS    SHOULD    FOLLOW    J08    CARD,  EACH  C 

C    BEGINNING    IN   COLUMN    1.  C 

c  c 

C         //     EXEC     FCRTCLGV,  REGION.  GO  =  350K  C 

C  / /FOR  T.SY  SPRINT    DD    DUMMY  C 

C  //SYSIN    DD    *  C 

r 

THE    SECONO     GF    THESE    SHOULD    5E    REMOVED    IF    A    PROGRAM  C 

C    LISTING    IS    DESIRED    AS    WELL    AS    OUTPUT    OATA.  C 

C  C 

C  C 

IMPLICIT   REAL*8    (A-H,0-Z) 

REAL* 3    XS  (25  ),Xi_(25  ),  XU(25  ),X(  2),XDOT<  2),DELTA(  31, 
*W(3)  ,  «N(3)  ,  BETA  (3  )  »THA CUT (3001  ,3)  ,X DAT A (3001, 5  )  , 
1A*°(3  ),0ELAY(3) 
D I  MENS  ICN    I  OUT  (3)  tlWTO  If  IT  AB (3 ) , INPUT ( 3  ) 
INTEGER     3TE°,RAMP  ,PA»A 

DATA    STE°,^  AM0t3ARA/i-HSTE0,  ^H°AMP  ,4HPARA/ 
COMMON    THAOUT,XDATA.T,0T,rF, 
*NQUT,  I  IN  ,NV  ,M3f  ICONT,  NEC, I  SKI  P,ITFfIOUT  ,NCPT 
COMMON     /REG1/    W, IPM 
C  C 

C  INPUT    CONTROL    FLAGS  C 

C  c 

REAO(  5,  SO  IN  RUNS,  N  OPT,  NGRAPH,NREF,  IIN,NOUT,  I  PRINT,  IFREQ 

WRITE  (  6,51  )NRUNS,  IIN,NOIJT 

IF<N0PT.EQ.1)WRIT  E(  6,  52  J 

I=(N0PT.EQ.2)  WRIT£(5,53) 

IF(NGRA?h.=0  .1  )WR  ITE(  6,  54) 

IF(NGRAPH.EQ,2)  WRIT  5(6  ,55) 

IF(NGRAPH.EQ. 3) ^R  ITE(  6, 56) 

IF(NREF.EQ.l)  WRIT  E(6,  57  ) 

IFMR£f=.EQ. 2)  WRITE(S, 58) 

IFdPPINT.EO.l    )WRITE(6,60)     IF^EQ 

IF(  IPRINT.EC2     )WRITE(6,59) 

^0    3  0    IPUN«1  ,NRUNS 

3EAD(  5,61  )T47,0T,  TF 

ITF=TF/DT 

IF(  ITF.GT.3001)  I  TF  =  3C0l 

iDP^ITF 

NPPLT=I0P/5 
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C 


■P^ 


WRIT5(6,62)T47,0T,TF 
C  C 

C  INPUT    OPTIMIZATION    DATA  C 

c  c 

GO    TO     (1,2),N0PT 

1  READ  (5,50 INV,NAV,NTA,NPR, IP,  IPM 
READ<  5,61  )<  XS(I  ), 1=1  ,NV) 
READ(5,61  )(XU(  I  ),I=I,NV) 
READ(5,61  XXL  (I  ),I=1,NV) 

WRI  TE( 6,6  3)NV,NAV,NTA,NPR,I P,I  P^ 
DO    31     1=1, NV 

WRIT  =  (6  ,64)  I,XS( I), I, XU (I ), I, XL (I) 

31  CONTINUE 

2  GO   TO    (3,4),NREF 

3  DO    32    I=1,N0UT 

C  C 

C  INPUT/CALCULATE    REFERENCE    CURVES  C 

C  C 

READ  (5,  50)  I  OUT  (  I)  ,  IWT{  I)  ,  I  TABU  ) 

W(  I  )sIWT( 11/1 00000. CO 

WRITE(6,65)     I,IOUT(I),I,     W     (I) 

IF(  IT*B( I) .EQ .1 ) WRITE ( 6,76) I 

ZF(irAB(I).E0«2  )WR  ITE<6,77  )I 

ITA=ITA3(I) 

GO    TO     (5,6),ITA 

5  RE AD (5, 61  )(XDATA( I  DAT  A, I) ,  IDATA=1,  ITF) 
GO    TO    32 

6  °EAD<5,75  )3ETA(  I  )  ,DELTA(  I) ,WN( 1 19 AMPC I ) yDELAY( I ) « 

*  I  N  P  UT  (  I  ) 

rfRITE(  6,66)  I  ,3ETA(  I  )  , I ,DE L^A ( I )  ,1  ,  «N  ( I  )  ,1  ,A«P(  I  )  , 

I,  INPUT  (I),  I,  OELAY(  I  I 

T  =  T47 

A17=0.D0 

IDATA=0 

X(l  )  =  0.D0 

X (2)=0.D0 

NT  =  0 

DO     33    J  alt  ITF 

11  IF(^.GE.D£LAY(  I  )  .AND.INOUT(  I)  .EQ.STEP  »A17=AMP<  I  ) 
ie(T.GS.DcLAY(I).ANO.  INPUT  t  I  )  .  EQ.  R  AMP  )  Al  7=AM°  (  I  ) 

*  *  ( T-0  EL  AY  {  I  )  ) 

IF(7  .GE.DELAY(  I  )  .AND.INPUT(  I  )  .EQ.PARA  )A17=AMP(  I  J 

*  *  (  (T-DELAY(  I)  )**2) 
XDQT( 1)=X(2) 

XDDT(2)=-2.C0*CELTM  I  )*WN(  I  )*X(  2  )-(WN  (I  1**2)  * 

*  X(1)+BETA<I )*A17 

S  =  RKLDEQ(  2,X,  XD  OT,  T,  DT,  NT  ) 

IF(S-1.D0  110, 11,16 
10  WRITE  (6,6  7)     I 

STOP 
16        XDATA<  JiI)*XU) 

33  CONTINUE 

32  CONTINUE 
GO    to    12 

4  DO     34    1=1, ITF 

DO    35    J=1,N0UT 
XDATA( I , J)=O.DO 
35  CONTINUE 

34  CONTINUE 

12  T=T47 
I3KI°  =  0 

C  C 

C  IP     SIMULATION    ONLY    CALL    PLANT  C 


C 


IF(  N0PT.EC2  )     CALL     3LANT(C) 
IF(N0PT.EQ.2)    GO    TO    13 


C 
C  OPTIMIZATION  C 

c  c 

5  =  1  .DO/3  .DO 
WRITE  (6,73  ) 
CALL  BOXPLX(NV,NAV,  NPR,NTA,R,XS  ,  I  =>,XU  ,  XL,  Ymnt  I ER) 
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WRITE(6,73) 
WRITE(6,63) 
DO    36    I'ltNV 

WRITEC6»69)If  XS(I  ) 

36  CONTINUE 

WRITE  (6,70)     YMN,IER 
C  C 

C  PRINTEC    OUTPUT  C 

C  C 

13  GO    TO    (I4115)f IPRINT 

14  WRITE<6,71) 

00    37    I  =  I,ITF,IFREQ 
TP=DT*I 
WRITE  (6, 72)  TPt(XOATMItJ)  ,THACUT(I,J)  ,J=1,N0UT  ) 

37  CONTINUE 

c  c 

C  PLOTTED   OUTPUT  C 

C  C 

15  GO    TO     (7,9,9)  ,NGRAPH 

7  03     38    1 1=1, NO  UT 

WP  ITE(6,73) 

CALL    o°L7(NPPLT  ,IDP,  II  ) 
33         CONTINUE 
GO    TO     9 

8  DO    39     11=1  ,NOUT 

WRITE(  6,73) 

CALL    PIC(NPPLTtlDPfll) 
39         CONTINUE 

9  WRITE(6,74)     IRUN 
30    CCNTINUE 

STOP 

50  FORMAT ( 815) 

51  FCRMATfO'  ,  12, IX,  •R'JN(5  )    WILL    *E    *AD6 ' , 4X, I 2» IX, 
* '  INPUTS' j4X,I2tlX,1  OUTPUTS'  ) 

52  FORMAT! «0', 'OPTIMIZATION   CALLED    FORM 

53  FCRMAT( '0*  , 'SIMULATION    CNLY    CALL6C    FOR') 

54  FORMAT! '0' , 'PRINTER     PLCT(S)     CALLEC    FOR') 

55  FCPMAT (  »0« ,  • VERS  AT  EC    PLOTS    CALLED    FOR1) 

56  =ORMAT{  '  0'  ,  •  MO    GRAPHICAL    OUTPUT    CALLED    CCR  •  ) 

57  FORM  AT  (  "0',  'REFERENCE    CURVE    CALLED    FORM 

53  FCRMAT < «0« , ' NO    REFERENCE    CURVE    CALLED    FOR') 

59  FGRMAT( '0' , 'NO    PRINTED     SIMULATION    DATA    CALLED    FOR') 

60  FORMAT(  «0S 'SIMULATION    DATA    PRINT    FREQUENC Y  =  '  , IX, I  2 ) 

61  FORMAT*  8  F10..5) 

62  FORMAT?  '0»  , 'INI  TIAL    TI  ME*'  ,F1  0.  5,  3X, '  STEP    SIZE*1, FID. 5 

*  ,  3X  ,  '  F  I  N'AL    T  IME-  •  ,  Fi  0  .5  1 

63  CCRMAT(  '0'  ,  'NV='  ,  I2,3X,'  NAV='  ,  12  ,3 X  , ' NT A= • ,  15 , 3X , ' NPR  = 
*»,I5,3X,  *I°  =  ,,Uf  3Xf  •  IPM*1  ,11) 

64  FCRMAT  (  '0' ,  'XS(  ' ,  12, •  ) =•  ,F1 0 .5, 5 X ,  •  XU(  » , I  2,  ' )  =  ' , F10. 5 , 
*5X,«XL(  •  ,12,'  )='  ,F10.5) 

65  FORM  AT  {  '0'  ,  '  IO'JT(  •  ,11,  •  )=• ,12, 5X,  'W(  ',!  1,' )  ='  ,F10. 5) 

66  FCRMAT  (  »0»  ,'  8ETA(  '  ,  II,  •  )=•  ,=10.5,  3X ,  '  DEL  TA  (  S  1 1 ,  •  )  =  •  , 
*F10.  5,3X,« WN(  ',11,'  )='  , FID. 5, 3X,'  AM  P(», II,'  )=•  ,F10  .5, 
*3X,  •  INPUT(  '  ,  II,  '  )  =  ',        A4    ,  3X,  'DELAY*  '  ,  I  1,'  )='  ,F10.5) 

67  cORMA~( '0'  ,'  INTEGRATION    TROUBLE-REFERENCE    CURVE') 

68  FORMAT*  'OS  'THE    FREE    PARAMETERS    ASSOCIATED    WITH    THE 
^MINIMUM    PERFORMANCE    MEASURE    ARE:') 

69  FCRMAT('Q»T'XS(',I2,')='  ,F20.10) 

70  FO'MiT (  '0? , 'THE    MINIMUM    PERFORMANCE    MEASURE    IS:',2X, 
*F20.10,5X,  '  IE?=«  ,2X,  II  ) 

71  FOR*AT(  «  1'  ,3X,'  TI  ME'  ,     9  X  ,'  XOA^A  ( 1  )  •  ,    7X  ,  '  TH  AQU~  (  i  )  '  , 

*  6X,  'XDATAt  2)  '  ,    7  X,  •  THACUTt  Z)  '  ,     6X  ,  '  XDA  ta(  3)  '  ,     "'X,  ' 
*ThA0UT(3  )'  ) 

72  FORMAT  (  '0'  t  6(  F10.  5,5X)  ,F1D.5) 

73  FORMA' CI') 

7«+    PGRMAT(  '  0'  ,' RUN    NUM8E  R  '  ,1X  ,  II  ) 

75  FORMAT (5  FID .5, A  4) 

76  FORMAT ('0'  ,«  REF.CURV E( ',  II, ')  =  ',  ITA«ULATED     DATA') 

77  FORMAT*  '0'  ,  'REF.    C  JRVE  (  ♦  ,11  ,  ■  )  =  S  EC  CND    ORCER    RESPONSE') 
END 
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SUBROUTINE  PLANT (C) 
C  C 

C         SUBROUTINE  PLANT (C)  SIMULATES  THE  SYSTEM  C 

c  c 

IMPLICIT    REAL* 3     (A-H,0-Z) 

REAL*3     G(25),R<  25  >,Z(  25), DRIVE  <25),THA<  25)  TOMG(  25)  , 
1THADOT(23),CMGDOT(25)  , CPV IN ( 25 ) , p LAG( 25 , 25  I , 
2    X2(  25)  ,    X2DOT(25),     THACUTO  001  ,3  )  , 
3XDATA(3001,3),C(25) t  AM P ( 25 ) . DEL A Y ( 2  5 1 , TAG ( 25,  25) 

DIMENSION    ICC251 , ID( 25  ),  IE (25 )  ,NF (25) ,IR (25)»  IV<25  )t 
1IEE<25),IMPUT( 25) ,I0UT(3) 
INTEGER    STEP,  RAM©,  PARA 

DATA    STEP,RAMP,PARA/<+hSTEP,4-HRAMP,4-HPARA/ 
COMMON    THAOUT,XDATA,T,DT,TF, 
CNCUT,  IIN,NV,M3,  I CONT  ,N EQ, ISK 1°, ITF,  I0UT,N0PT 
IF    (ISKIP-1)    1,5,5 
1    I  SKIP    =    2 
REA0(5,24>    N 
HI*    OT 
h2   =    0.5D0*Hl 
Nil     =    0 
N55    =     0 
M66    -    0 
ICK    =    2*N 
WRITE  (6,310) 
DO    3    1  =  1, N 
C  r 

C  INPUT    BLOCK   OATA  C 

C  C 

READ (5,25  )  IC ( I ) ,  ID (  I) ,  IE (  I  ) ,  IV (  I)  ,G( I  ) , P ( I ) , Z (  I ) 
IF<  ID(  I  ).  EQ.l)     Nil  =  NllH 
IF(  ID(  I  ).EQ.5  )    N55=N55+1 
IF(  IQ(  n  .  EQ.6J     N66-N66*! 

WRITE (5, 26)     IC( I) ,10(1 ), IE(  I ) , IV( I),G(  I ) ,P(  I ) ,Z(  I  ) 
3         C  CNT  I  NU  E 

WRITE  (6,312) 
00    316    1  =  1,  UN 

c  c 

C  INPUT    FORCING    "UNCTICNS  C 

C  C 

READ( 5,311 )AMP{ I),DELAY( I), INPUT!  I),  IEE(  I  ) 
WRI TEC  6,317)  I SE(  I )  ,AMP(I)  , INPUT (I) , DELAY (  I  ) 
316  CONTINUE 

Nil     =    4* Nil 

N66    =    4*N66 
NEQ    =    N-l 

C  C 

C  SET    POINTERS  C 

DO    4    J=1,N 
DC    4    K  =  l, N 

CLAG  (K,  J)  =O.D0 

IF(  I  V(K).EQ.  IE(  J)  )F  LAG  (K,  J)  =1  .DO 
<*         CCNTINUE 
C  C 

C  INITIALIZATION  C 

C  C 

5    DO    6    ICLR=!  ,N 

THAI  ICLR  )  =  0.D0 
THADOT( ICLR )  =0  .  00 
OMG( ICLR)=0. DO 
OMGOQT (  ICL'  )=0.D0 
X2  (  ICLR)=0.  DO 
X2D0K  ICLR)  =0.D  Q 
6  CCNTINUE 

T=O.CO 
NR2    =    1 
NR3    =    1 
NR4    =    1 
M  3   =    0 
I CONT    =    0 
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IWAIT     =    0 

IT   =    0 

ILAST    =     0 

I5LAST    =    0 

I6L4ST    =    0 

IF(^3PT  .EQ  .2)    GO    TO    7 
C  C 

C  ENTER    C(I)=G(I)     VALUES    AT    THIS    PCINT  C 

c  c 

7    DG    9    MDRV=1,N 

DRIVE(MDRV)=O.DO 
DRVIN  (MDRV)=O.DO 
00    8     M=l  ,N 

IFt IV<M)„NE.IE(MDRV))    GO    TO    315 

DRIVE(MCRV)  =  DRIVE(MCRV  )  +T  HA  (  M  )  *FLAG  (  M,  MOR  V) 
3  15  I<=(IEE(M).NE.  IE(MOPV))    GO    TO    8 

I«=  (  INPUT(M)  .  EQ.STEP)    GO    TO    32 

IF( INPUT (Ml .EQ. RAMP)    GOTO    32 

Ic ( INPUT(M) .EQ.PARA)    GO    TO    34 

32  IF(T.GE.OELAY(M ) ) DRVIN (MOR VI =AMPtM) 
GO    TO    9 

33  IF(T.GE.DELAYM))  ORVIN(MORV  )  =4MP(  M  )  *  (  T-DELAY  (M)  ) 
GO   tp    3 

34  I F (T.GE. DELAY (M) ) DRVIN (MDRV )»  AMPtM  )  *UT-OEL  AY(  M  )) 
*           **2) 

3  CONTINUE 

OR  I  VE(MORV)  =ORIVE(  *DRV  )  +ORV  I  N  (M[RV) 
9  CONTINUE 

10  M3    =    M3+1 

:  c 

C  BLOCK     SIMULATION  C 

c  c 

IP  ( IWAIT.SQ.OJ  T  =  T+H2 
IF  (I0(«3)  .SQ.ll  GO  TO  11 
IF  (  I0(  M3  )  .E0.2  )  GO  TO  12 
IP  ( ID(M3J.=0.3)  GO  TC  13 
IF  (  ID( ^3  ).E0.4)  GO  TO  14 
IF  (I0(  M3)  .EQ.5  )  GO  TO  15 
IF  ( 10 ( M3) .EQ.6)  GG  TO  16 
WRITE  (6,23) 
STOP 

11  THAM31     -    G(M3)*DRIVE(M2) 
IWAIT     =     IW4lT-t-l 

ILAST    ~    ILAST-t-1 

IF    U  ILAST. EQ. Nil)  .ANC.  (ICONT. EG. NEC)  )    GO    TO    21 

r  r     t  n      1  Q 

12  THA0QT(M3)     =    -P (M3)*THA (M3 I +G( M3 ) *CPIVE (M3  ) 
S    =  RKLDE2<  THA,  THA00TtNR2) 

IWAIT    =    IWAIT+1 

I  F(  S-1.DGJ17, 13,21 

13  0MG00T(M3)     *   -P(M3)*0MG(M3)+G(M3)*0RIVE<M3) 
S     =    RKL0E3  (  CMG.  OMGOOT  ,NR3) 

THAM2)     =    (Z(M3)-P(M3n*GMG(M3  )  +  G(M3)*0RIVE(M3) 

IWAIT    =     IWiIT  +  1 

IF(S-1.00)17,18,21 

14  V    -   CC PLX< °,Z,G, DRIVE, X2,0MG,NR4) 
THA( M3)     =    0MG(M3 ) 

IWAIT    =    IWAIT+1 

I  F{  7-1.  DC)  17,  13,21 

15  THA(*3)     =    G(M3)*ORIVE(M31 

Ic  (THA(M3)  .GT.P(M3)  )  ThMM3)=P(*2) 
IF  (THA(M3)  .LT  ,Z(  M3)  )  THA  (  M3  )  =  Z  (  M  2  ) 
IWAIT     =     IWAITH 

I  3LAST    =    I5LAST+1 

IF    (  (  I5LAST.EQ.N55)  .AND.  (ICONT.EQ  .NEO)  )     GO     T0    21 

GO    TO    18 

16  THMM3)     =    G(M3)*0RI  VE(  ^2) 

I  F((  DaBS  (T  HA  (  M3  )  )  )  .LT.PM3)  )  THA  (  M  I )  =0  .0  0 

IWAIT    =    [WAIT+1 

loLAST     =    I6LAST*! 

IF    l(I6LAST.EQ.N66).ANC.(  ICONT  .EQ  .NEQH    GO    TO    21 

GO    TO    1 3 
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17    WRITE     (6,29) 

STOP 
13    ICONT    =     ICDNT+l 

IF    (ICONT-N)    10,19,20 

19  MR2    =    NR2+1 
NR3    =    NR3+1 
NR4    =    NR4+1 
M2   =    0 
ICONT    -    0 

IF    (IWAIT.EO.ICO     IWAIT»0 
GO    TO    7 

20  WRITE    (6,30) 
STOP 

C  C 

C  STORE    OUTPUT    DATA  C 

C  C 

21  IT  »    IT+1 

OC   33    I»ltNCUT 

THAOUT(IT,I)  =THA(IOUT(I)  )       ' 
38    CONTINUE 

IF    (T-TF)     22,22,23 

22  MR2  =  1 
NR3  =  1 
NR4    =    1 

ICONT    =     0 
I  WAIT    =    0 
ILAST    =     0 
I5LAST     -    0 
I  6LAST    =    0 
GO    TO    7 

23  RETURN 

24  FORMAT*  12) 

25  FORMAT  (412,  2X,3F10.5  ) 

26  =  CRMAT     (/,5H    BLK    ,I2,I1,2H    =  ,2  12  ,  6X,  3  520  ,7  ) 

27  FORMAT     (//, 2X,'TH£TA    OUT    IS    TH  A  (  •  ,  I  2  , «  )  '  ) 

23    FCRMAT     (40H    ***    EQN    SWITCH    CONTROL    DID    NOT    WORK    ***) 

29  FORMAT     (20H    INTEGRATION    TRQU8LE) 

30  F0RV»AT(58X,  'ERROR     IN     INTER  GATI  ON.  ',  2X  , 
1' ATTEMPTED    TO    INTEG.     MCRS    THAN    N-EC^N*  ) 

310  FORMAT ( io«f40X,» THE     SYSTEM    BLOCKS    ARE:') 

311  FCPMAT  (2F10.5,A4,  12  ) 

312  FORMAT (' C  ,40X,  •  THE    SYSTEM    FORCING    FUNCTIONS     ARE:') 
317    FORMA"  (  »0',  'DRV  IN  (  «,  12,  '  )=  '  ,  F  10.  5  ,  A4,F  10.  5) 

END 

FUNCTION    RKL0£Q(N,X,X0OT,T,0T,NT) 

C  C 

C  CALCULATES    RESPONSE    OF     SECOND  C 

C  ORDER    REFERENCE    CURVE  C 

c  c 

IMPLICIT    REAL*8    (A-H,0-Z) 
REAL*8    X(2)  ,XDOT(2),0(25  ) 
NT  s    NT   ♦  1 
GO    TO     (1,2,3,4) ,NT 

1  HI   *   DT 

H2    =    H1*C.5D0 
H3    =    HI  *2.0D0 
H6    -    H1/6.0DO 
DC    11    J    =    1,N 
11    0(J)     =    O.DO 
A    =    0.5D0 
T    =    T     *•    H2 
GO    TO    5 

2  A    =    0,2923932133134525 
GO    TO    5 

3  A    =1.7071067311365475 
T    -    T    «■    H2 

GO   TO     5 

4  DC    41     I    =    1,N 

41    X(I)     -     X(I)     4-    H6*XD0T(I)    -    Q(I)/3.D0 
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NT   »     0 

RKLDEQ    =    2. 
GC   TO    6 

5  00    51    L   =  1,N 

X(L)    =    X(L)     ♦    \*i  0T*X0QTtL)-0U)  ) 
51    Q(L)    =    H3*A*XD0T(L)     *     ( 1 .  CO    -    3  •  00*A  )  *Q  (  L  ) 
RKLQEQ     =    1. 

6  RETURN 
END 


FUNCTION    RKLDE2     (X,XD0T,NR2) 
IMPLICIT   REAL*3    (A-H,Q-Z) 
REAL*3    THAOUT(3001,3)  , XOATA (  3001,3) 
REAL*3    X(4),    XD0T(4),     G(25) 
DIMENSION    I0UT(3) 
COMMON    THA0UT,XOATA,T,DT,TF, 
CNCUT,II  N,NV,M3,ICaNT,NEC,ISKI=>,  IT  F,  IOUT 
GO    TO     (  1,2,3,4)  ,    NR2 

1  Hi    =    OT 

H2    =    Hl*0.500 
H3    =    Hl*2.0D0 
H6    =    H1/6.0D0 
Q(M3)     =0.00 
A    =    0.  50  0 
GC    TQ    5 

2  A    -    0.2923932138134525 
GO    TO    5 

3  A    =    1.7071067811365475 
GO    TO    5 

4  X(M3)    =    X(M3)*H6*XOOT(M3)-Q(M3) /3.DC 
RKL0E2     *    1. 

IF    (  ICONT.EQ.NEQ)     RKLDE2=2. 
GO    TO    6 

5  X(M3)     =    X(M3)+A*(  CT*XDCT(M3)-Q(M3  )) 
Q(M3)     -    H3*A*XD0T<  M3)-M1.00-3.DO*A)*Q<  M3) 
PKL0E2    =    1. 

6  RETURN 
ENO 


FLNCTION    RKLDE3     (X,XDCT,NR3) 

IMPLIC  IT    0EAL*6     (A-H,0-2) 

REAL* 3     THAOUT(3  001,3  )  , XOATA (3001 ,  3) 

REAL*3    X(4),    X0QT(4)f    Q(25) 

0IM6NS  ICN    I  CUT  (3) 

CCMMCN     THAOUT,  XOATA,  T,CT,T=, 
CNOUT,  IIN,NV,M3,  ICONT,NEQ,I SKIP , IT F,  IOUT 

GC   TO     (1,2,3,4),    NR3 
,    HI    =    OT 

H2    -     H1*0.5D0 

H3    '    Hi  *2.0  00 

H6    =    HI  /6.  000 

Q(M3 )    =    0.00 

A    =    0.5  00 

GO   TO    5 
:    A   =    0.2923932133134525 

GC   TO    5 

A    =    1.7071067311365475 

GC   TO    5 
•    X(M3)     =    X(M3)+H6*XDOT(M3)-Q(M3)/3«D0 

RKLDE3    =    1. 

IF    (  ICQNT.EC.NEQ)     RKL0E3=2. 

GC    TO    6 
i    X(*3)     =     XM3)+A*<  0T*XD0T(M3)-Q(  M2)  ) 

Q(M3)     =    H3*A*X00T  (M3)+<i.00-3  .00*A)*0(M3J 

RKL0E3    =    1. 
i    RETURN 

ENO 


FUNCTION    CCPLX     (P  ,  Z,  G,  CR  IV  E  ,  X2  ,  OMG,  MP4  ) 
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IMPLICIT    REAL*8     (  A-H,  O-ZI 

REAL*  3    THAOUT  (3001,31  »XCATA(  3  00'  ,3  }  r CMGDGT  (1) 

REAL"* 8    0MG(l)«P<l)f  Z(l)tG(l)fORlVE(I)  fX2(25TtX20C7(25) 

01  MENS  I  ON    I0U.T(3) 

COMMON    THAOUT,XDATA,T,OT,TF, 
CNOUT,  IIN,MV,M3,  ICONT,NEQ,  I  SKIP,  ITF,IOUT 

QMGD0T(M3J     =   X2 (M3) 

SS    =    RKLDE3(0MG,QMG00T,NR4) 

X2D0T(M3)    =    -2.*P(M3)  *Z(M3)*X2(  M2)-Z(M3  )**2*QMG(M3)  + 
1G(M3)*DRIVE(M3 ) 

SSS   =    RKLDE4(X2,X2D0T,NR4) 

CCPLX    =    SSS 

RETURN 

END 

FLNCTION    RKLDE4    <X,XDCT,NR4) 
IMPLICIT    REAL*8    (A-H, 0-Z) 
REAL* 8    THAOUT  (3  001,3  ) ,  XOAT  A(  3001 ,  2  ) 
REAL*8    X(l)  ,    XOGT(l),     CC(25) 
DIMEMS  ION    IOUT( 3) 
CC*MCN    THAQUT,X  DATA  ,T  ,  DT ,TF, 
CNOUT, II N,NV,M3,IC0NT,NEQ,ISKIP,ITF,I0UT 
GO    TO     (1,2,3,4.)  ,    NR4 

1  HI    =    OT 

H2  =  H1*Q.500 
H3  =  HI  *2.0D0 
H6  =  Hl/6.000 
0C(M3)  -  0.00 
A  =  0 .500 
GC   TO    5 

2  A    =    0.2928922188134525 
GC   TO    5 

3  A    =    1.7071067311365475 
GO   TO    5 

4  X(M3)     =    X(  M3  )-H-i6*X0OT(M3)-QC(M3)/2.0O 
RKL0E4    '    1. 

IF    (ICONT.EQ.NeQ)    RKLDE4=:2. 
GC   TO    6 

5  X(M3)     =    X(M3)-»-A*(0T*XQQT(M3)-qC(  «3 )  ) 
QC(M3)     =    H3*A*XO0T(  M3  )*<  1.00-3. DO*A)*QC(  M2) 
RKL0E4    =   1  . 

6  RETURN 
ENO 

FUNCTION    KE<C) 

c  c 

C  EVALUATES    IMPLICIT    CONSTRAINTS    FCR    3CXPLX  C 

C  C 

REAL*3     C(25) 

KE=0 

R  ETURN 
ENO 

FUNCTION    FE(C  ) 

c  c 

C  THE    FUNCTION    MINIMIZED    BY     SCXPIX  C 

C 
IMPLICIT    REAL *3  (A-H, 0-Z  ) 
REAL* 3    THAOUT (3 00 1,3)  ,XDATA(3 001  ,3) f C (25 ) , W (3 ) 

COMMON     THA0UT,X0A7A,T ,OT,^F, 
*NOUT,I  INtNV  ,M3,ICQNT  ,  N5CISKIP,  ITC,  InUT 
COMMON     /PEG1/    WfIPM 
DC    ID     I=1,NCUT 

THAOUT  (1  ,1)  =0.0  0 
10    CONTINUE 
OIFF=O.CO 
PI=0.00 
CALL    °LAN"(C) 
T=0T 
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DC   il     1  =  1, ITF 
DO    12    J=l  ,NOUT 

3IFF=XD4TA(  I,  J  )-THAQUT(  I  ,  J  ) 
GO    TOd  ,2,3,4-  ),  IPM 

1  PI  =  PI*W  (J)*(DIFF*2) 
GO    TO    12 

2  PI=PI*W  (J)*(DIFF**2  J*T 
GO    TO    12 

3  °I  =  PI+W  (J)*DA3S(D  IFF  ) 
GO    TO   12 

<►  PI  =PI  +  W(  J)*0A3S(DIFF  )*T 

12         CONTINUE 

T=T+DT 
11    CONTINUE 
FE==>  I 
RETURN 
END 

SUBROUTINE    PPLT  (NPPLT  ,  IDP,  1 1  ) 
C  C 

C  AUTOSCALES    AND    CALLS    FOR    PRINTR    °L3TS  C 

C  C 

IMPLICIT   PEAL*8     (A-H,C-Z) 

P.EAL*4-    XX( 900), YY<900) ,WW(  900) 

REAL*4.       TX(  4)  ,TY(  4-) 

REAL*  3    THAOUTOOOlt  3)  ,  X0AT4 ( 3 001 ,3 ) 

DIMENSION    ICUT(3) 

CCMMGN    THA0UT,X3ATA,T,DT,TF, 
CNOUT,  IIN,NV  ,M3  ,ICONT,  NEC,  I  SKI  P,  ITF,  I  OUT 

DC    1     1=1,900 

XX(  I  )     =    0. 

YY( I )     =    0, 

1  WW< I )     =    0. 
T  =  0.  DO 

TSTEP  =  5.0Q*0T 

J    =    0 

3IGX=Q.D0 

3  I GY  =  0.00 

SMLX=O.DO 

SMLY=O.DQ 

DC    2     1  =  1, IDP, 5 

J    =    J  +  l 

XXf  J  )     *   T 

YY< J)=XOATA(  I, II  ) 

WW<  J)  =THAOUT(I  ,11  I 

X    =    T 

XD   =    YY (J) 

TH    =    WW< J) 

YMAX=DMAX1(XD, TH> 

YMIN  =  D.MIM  (XOtTHI 

XVAX  =  DMAX1 (3IGX  ,X) 

I  F    (3  IGX.LT.X)    3IGX=X 

IF    (3IGY.L~.YMAX)     BI3Y=Y*4X 

IF    (  SMLY.GT.YMIN)     SMLY*YMIN 

T    =    T+TSTEP 

2  CCNTINUE 
TX(1)     =    0. 
7X(2  )    =    0. 
TX(3)     =    SMLX 
TX(^)     =     3IGX 
TY(1 )     =    3I3Y 
TY(2)     =    SMLY 
TY(3)    =    0. 
TYm     =    0. 

WRITE     (6,3)     3IGX, 3IGY.SMLY 

C&LL    °LOTo      (TXtTYf4.tl  I 

CALL    PLC'°     (XX, WW, NPPLT  ,2) 

CALL    OLOTP     ( XX, YY,NPPLT,3) 

WRITE  (6.4) 

I P( IDP.EQ.4500)    W.RITE(6,5) 

'ETUPM 


3  FORMAT* 2X, • BIGX=    •  ,  El  5  .  7  ,2  X,  •  B  I  GY  =    ',E15#7»2X, 
l'SMLY=     *,E15.7) 

4  FORMAT* //,2X,» SYSTEM    RESPONSE    FOR    PRCBLEV    ') 

5  FORMAT  ( //,2Xt.'STDP    AT    900    GRAPH    POINTS') 
END 

SUBROUTINE    PIC    (NPPLT ,  IDP,  1 1  ) 
C  C 

C  AUTOSCALES    ANO    CALLS    FOR    VERSATEC    PLOTS  C 

c  c 

IMPLICIT   REALMS    (A-H,C-Z) 
PEAL*<V    XX(  900),YY(90Q)  ,V»W(  900) 
REAL*4TX(4)  ,TY(  4) 

REAL*3    THA0UTO001,  3)  ,  XOATA ( 3 001 , 2 ) t TI  TLE (12 ) 
REAL   *3LABC/'  •  / 

REAL    *4LAB4/'     A       »/,LABC/'        0    «/ 
DIMENSION    I0UT(3) 
CCMMON    THAOUT,X0ATA,T  ,DT,TF,     ' 
*MQUTtIIN,NV  ,M3»ICaNT,  NEC  ,1  SKI  => ,  IT  F,  I  GUT 
REAO    (5,2)    TITLE 
TO. 00 

TSTEP=5.J0*0T 
J    =    0 

3IGX    =    0.00 
3IGY    =    0.00 
SMLX    =    O.DO 
SPLY     =    0.00 
00    1    I«l,IOPfS 
J    =    J+l 
XX(J)      =    T 
YY(  J)  =  X0ATA  (1,1  I  ) 
WW  (J  )=THAQUT(  I,  II  ) 
X    =    T 

XO    =     YY( J) 
Th    =    WW (J) 
YMAX    =     CMAXl(XO,Th) 
YMIN    '    0MIN1 (XO,TH) 
XMAX    =     CM  A  XI  (  3IGX,X  ) 
IF    (3IGX.LT. X)     3IGX=X 
IF    (BIGY.LT.YMAX)    8IGY*YMAX 
IF    (SMLY.GT.YMIN)    SML*=YMIN 
T    =    T+TSTE3 
i    CONTINUE 
TX(1)     =    0. 
TX(Z)     =    0. 
^X(3)     =    SMLX 

Txm    =   3iGx 

TY(1)     =    BIGY 

TY(2)     =     SMLY 

TY(3 )     =    0. 

TY<4)     -    0. 

CALL    DRAW    (4,TX,TY,1, lf LABCt TITLE t C, 0, 0, C,  C,  0,  3  1 3  ,  0  ,L  ) 

CALL  DRAW  (NPPLTtXXfWW,2,0»LABAtTITLE,OtO,OfO,0,0,8,8, 
*0,L) 

CALL  DRAW  (NPPLT, XX, YY,3tO,LA3D,TITLE,  G,C,0,0,O,O,3,3, 
*0,L) 

IF    (L.NE.O)     WRITE     (5,2)    L 

RETURN 

2  FORMAT     ( 5A3) 

3  FORMAT     (//,•       GRAPH    NOT    COMPLETED.    OUTPUT    COOE     =    ',12) 
END 

SUBROUTINE    30X  =  LX     ( NV , NAV , NPR , NT Z , R Z, XS ,  I P , 3U , 3L , YMN  , 
*IER) 

c  c 

C  SUBROUTINE    30XPLX  (CATEGORY    HO)  C 

C  C 

C  PURPOSE  C 

c  c 
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C  30XPLX    IS    A    SUBROUTINE   USED    TQ    SOLVE    THE     PROBLEM  C 

C  OF    LOCATING    A    MINIMUM     (JR    MAXIMUM).  OF    AN     ARBITRARY  C 

C  OBJECTIVE    FUNCTION    SU3JECT    TO    ARBITRARY    EXPLICIT  C 

C  AND/OR    IMPLICIT     CONSTRAINTS     BY    THE    COMPLEX    METHOD  C 

C  OF    m.j.     BCX.     EXPLICIT    CONSTRAINTS    ARE    CEFINED    AS  C 

C  UPPER    AND    LOWER    30UNCS    ON    THE     INDEPENDENT    VARIABL-  C 

C  =  S.     IMPLICIT    CONSTRAINTS    MAY     SE    ARBITRARY    FUNCTIQNSC 

C  OF    THE    VARIABLES.    TWO    FUNCTION    SUBPROGRAMS    TC  C 

C  EVALUATE    THE    OBJECTIVE    FUNCTION    AMD    IMPLICIT    CON-  C 

C  STRAINTS,    RESPECTIVELY,    MUST    BE    SUPPLIED    BY    THE  C 

C  USER     (SEE    EXAMPLE    BELOW).    30XPLX    ALSO    HAS    THE    OPT-  C 

C  ICN    TG    PERFORM    INTEGER    PROGRAMMING,    WHERE    THE    VAL-  C 

C  UES    OF    THE    INDEPENDENT    VARIABLES    ARE    RESTRICTED   TC  C 

C  INTEGERS.  C 

C  C 

C  USAGE  C 

C  C 

C  CALL    BCXPLX     (NV ,NAV ,NPR ,NTA, R , XS,  IP, XU, XL, YMN ,  IER  )  C 

C  C 

C  DESCRIPTION    OF    PARAMETERS  C 

C  C 

C  NV          AN    INTEGER     INPUT    DEFINING     THE    NUMBER    OF    INDE-  C 

C  PENDENT    VARIABLES    OF    THE    OBJECTIVE    FUNCTION  C 

C  T0    3E    MINIMIZED.     NOTE:     MAXIMUM    NV    +    NAV    IS  C 

C  PRESENTLY    50.    MAXIMUM    NV     IS    25.     IF    THESE  C 

C  LIMITS    MUST    Be    EXCEEDED,     °UNCH    A    SOURCE    DECK  C 

C  IN    THE    USUAL    MANNER,    AND    CHANGE    THE    CI  MEM-  C 

C  SION    STATEMENTS.  C 

C  C 

C  NAV       AN    INTEGER     INPUT    DEFINING    THE    NUMBER    OF    AUXI-  C 

C  LIARY    VARIABLES    ThE    USER    WISHES    TO    OEFINE    FOR  C 

C  HIS    OWN    CONVENIENCE.    TYPICALLY    HE    MAY    WISH  C 

C  TO    DEFINE    T riE    VALUE    OF    EACH    IMPLICIT    CON-  C 

C  STRAINT    FUNCTION    AS     AN     AUXILIARY    VARIABLE.    IF  C 

C  THIS    IS    DONE,     THE    CPTI ONAL    CUTPUT    FEATURE  C 

C  OF    BOXPLX    CAN     BE    USED    TO    OBSERVE    THE    VALUES  C 

C  OF    THOSE    CONSTRAINTS    AS    THE    SOLUTION    =>R0-  C 

C  GRESSES.    AUXILIARY    VARIABLES,    IF     LSED,     SHOULD  C 

C  BE    EVALUATED     IN    FUNCTION    K6    (DEFINED    3EL0W).  C 

C  NAV    MAY     3c     IERQ.  C 

c  c 

C  NPR       INPUT    INTEGER    CONTROLLING    THE    FREQUENCY    OF  C 

C  OUTPUT    DESIRED    FCP    DIAGNOSTIC    PURPOSES.     1=  C 

C  NPR    .LE.     0,    NO    OUTPUT    WILL    BE    PRODUCED    3Y  C 

C  BOXPLX.    OTHERWISE,     THE    CURRENT    COMPLEX    OF  C 

C  K  =  2*NV    VERTICES    AND    THEIR     CS\"~OI0    WILL    BE  C 

C  OUTPUT    AFTER    EACH    N=>R     PERMISSIBLE    TRIALS.    THE  C 

C  NUM3EP    OF    TOTAL    TRIALS,     NUMBER    OF    FEASIBLE  C 

C  TRIALS,    NUMBER    OF    FUNCTION    EVALUATICNS     AND  C 

C  NUM8ER    OF     IMPLICIT    COMSTRAINT    EVALUATIONS    ARE  C 

C  INCLUOED    IN    THE    OUTPUT.  C 

C  ADDITIONALLY,     (WHEN    NPR    ,GT.    0)     THE    SAVE    IN-  C 

C  FORMATION    WILL    BE    OUTPUT:  C 
C 

C  1)     IF    THE     INITIAL    POINT    IS   NOT    FEASIBLE.  C 

C  2)    AFTE^    THE    FIRST    COMPLETE   COMPLEX    IS  C 

C  GENERATED.  C 

C  3  )     IF    A    FEAS  I3LE    VERTEX    CANNOT    BE    FOUND    AT  C 

C  SOME    TRIAL.  C 

C  4)     IF     THE    OBJECTIVE    VALUE    CF    A     VERTEX    CAN-  C 

C  NOT    3E    MAGE    NO-LON  GER- WPP  ST  .  C 

C  5)     IF    THE    LIMIT    ON    TRIALS     (NTA)     IS    REACHED  C 

C  AND,  C 

C  6)     WHEN    THE    OBJECTIVE     FUNCTION    HAS     BEEN    UN-  C 

C  CHANGED    FOR    2*NV    TPIALS,     INDICATING     A  C 

C  L3CAL    MINIMUM    HAS     BEEN    FOUND.  C 

C  C 

C  Ic    THE    USER    WISHES    TO    TRACE    THE     PRCGRESS     OF  C 

C  A    SOLUTION,     A     CHOICE    OF    MPcs     25,      5C    OR    100    IS  C 

C  RECCMMENOSD.  C 

C  C 

C  NTA        INTEGER     INPUT    CF    LIMIT    ON    THE    NUMBER    OF  C 
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C  TRIALS    ALLOWED     IN    THE    CALCULATION.       IF     THE  C 

C  USER    INPUTS    NT  A    .LE.    0,     A    CEFAULT    VALUE    OF  C 

C  2000    IS    USEO.       WHEN    THIS    LIMIT    IS    REACHED  C 

C  CONTROL    RETURNS    TO    THE    CALLING    PROGRAM    WITH  C 

C  THE     BEST    ATTAINEC    OBJECTIVE    FUNCTION    VALUE    IN  C 

C  YMN,     AND    THE    BEST      ATTAINED    SOLUTION    °CINT    IN  C 

C  XS.  C 

c  c 

C  R             A   REAL    NUMBER     INPUT    TO    DEFINE    THE    FIRST    RAN-  C 

C  DGM    NUM8ER    USED    IN    DEVELOPING   THE    INI-  C 

C  TIAL   COMPLEX    OF    2*NV    VERTICES.  C 

C  (0.    GT.    R     .LT.     1.  )     IF    R     IS    NOT    WITHIN    THESE  C 

C  BOUNDS*     IT    WILL    BE    REPLACEC    BY    1./3.    .  C 

C  C 

C  XS          INPUT    REAL    ARRAY    DIMENSIONED    AT    LEAST    NV+NAV.  C 

C  THE    FIRST    NV    MUST    CONTAIN    A    FEASIBLE    ORIGIN  C 

C  FOR    STARTING    THE    CALCULATICN.       THE    LAST    NAV  C 

C  NEED    MOT    BE    INITIALIZED.       UPON    RETURN    FROM  C 

C  BOXPLX,     THE    FIRST    NV    ELEMENTS    OF    THE    ARRAY  C 

C  CONTAIN    THE    COORDINATES    OF     THE    MINIMUM    OB-  C 

C  JECTIVE    FUNCTION,     AND    THE    REGAINING    NAV  C 

C  (NAV    .GE.     0)     CCNTAIN    THE    VALUES    OF    THE  C 

C  CORRESPONDING    AUXILIARY    VARIABLES.  C 

C  C 

C  IP          INTEGER    INPUT    FOR    OPTIONAL    INTEGER     PRO-  C 

C  GRIMING.       IF     IP=1,     THE     VALUES    OF    THE    INDE-  C 

C  PENDENT    VARIA8LES    WILL     BE    REPLACED    WITH  C 

C  INTEGER    VALUES     (STILL    STCRED    AS    REAL*4).  C 

C  C 

C  XU          A    REAL    ARRAY    DIMENSIONED    AT    LEAST    NV     INPUT-  C 

C  TING    THE    UPPER     BOUND    ON    EACH    INDEPENDENT  C 

C  VARIABLE,     (EACh    EXPLICIT    CONSTRAINT.       INPUT  % 

C  VALUES    ARE    SLIGHTLY    ALTEREC    BY    BOXPLX.  C 

C  C 

C  XL          A    REAL    ARRAY    DIMENSIONED    AT    LEAST    NV    IN^UT-  C 

C  TING    THE     LOWER    8CUND    ON    EACH    INDEPENDENT  C 

C  VARIABLE,     (EACH    EXPLICIT    CONSTRAINT).       NOTE:  C 

C  FOR    EOTH    XU    ANO    XL    CHOOSE    REASONABLE    VALUES  C 

C  IF    NONE    ARE    GIVEN,    NOT    VALIES    WHICH    ARE    MAG-  C 

C  NITUDES    ABOVE    OR    BELOW    THE    EXPECTED     SOLUTION.  C 

C  INPUT    VALUES     ARE    SLIGHTLY     ALTEREC    SY    BOXPLX.  C 

c  c 

C  YMN       THIS    OUTPUT    IS    THE    VALUE     (?EAL*4)     OF    tH£    03-  C 

C  JECTIVE    FUNCTION,    CORRESPONDING    TO    THE    SOLU-  C 

C  TION    POINT    OUTPUT    IN    XS.  C 

C  C 

C  IER       INTEGER    ERROR    RETURN.       TO     EE    INTERROGATED  C 

C  UPON    RETURN    FROM    BOXPLX.        IER    WILL    3E    ONE    OF  C 

C  THE    FOLLOWING:  C 

:  c 

C  =-1      CANNOT    FIND    FEASIBLE    VERTEX    CR    FEAS-  C 

C  I  BL  E    CENTROID    AT    THE    START    OR    A    RE-  C 

C  START    (SEE    '"ETHOC     BELOW).  C 

C  =0         FUNCTION    VALUE    UNCHANGED    FOR     »N'  C 

C  TRIALS.        (WHERE    N*6*NV+10)       THIS    IS  C 

C  THE     NORMAL    RETURN    PARAMETER.  C 

C  =1         CANNOT    OE^ELO0    FEASI3LE    VERTEX. 

C  -Z         CANNCH"     DEVELOP    A    NO-LCNGER-WGRST    V  ER-  C 

C  TEX. 

C  =3         LIMIT    -)N    TRIALS    REACHED.        (  NTA     E  X- 

C  CEEDED) 

C  NOTE:        VALID     RESULTS    MAY    BE    RETURNED     IN  C 

C  ANY    OF    THE    ABOVE    CASES.  C 

C  C 

C  EXAMPLE    OF    USAGE  C 

C  C 

C  THIS    EXAMPLE     MINIMIZES    THE    OBJECTIVE    FUNCTION  C 

C  SHOWN     IN    THE    EXTERNAL    FUNCTION     FE{  X)  .        THERE    ARE  C 

C  TWO    INDEPENDENT    VARIABLES    X  (1  )     S    X(2),     AND    TWO     I M-  C 

C  PLICIT   CONSTRAINT    FUNCTIONS     X(3)     S    X  ( 4 )     WHICH    ARE  C 

C  EVALUATED    AS    AUXILIARY    VARIABLES     (SEE     EXTERNAL  C 

C  ^UNCTION    KE(X)     )  .  C 
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DIMENSION       XS  (4)  ,  XU  (  2)  ,XL(  Z) 

STARTING    GUESS    • 

XS(1)     =    1,0 

XS(2)  =  0.5 
UPPER    LIMITS 

XU(1)    =     6.0 

XU(2  )  ■  6.  0 
LOWER    LIMITS 

XL(1  )    =    0.0 

XL(2)     -    0.0 

R    =   9./13. 

NTA    =    5000 
NPR    =     50 
NAV    =    2 
NV   =    2 

1°    =    0 

CALL    60XPLX    ( N V, NA V  ,NPR  ,NTA , R , XS , I P, XU ,XL t YMN, I ER ) 
WRITE (6,1)     UXS(  I  ),  1=1,  4)  ,YMN,  I  ER) 
1    FORMAT*////  tf  THE     PCINT    IS    LOCATED    AT     (XS (  I »  = )    '  , 

*4(E13.7,5X) 9//t  •         AND    THE    FUNCTICN    VALUE    IS    •- 


*E13 
STOP 
END 


7t  • 


IER    = 


I  5) 


FUNCTION    KE  (X; 

EVALUATE  CONSTRAINTS.  SET  KE-0 
CONSTRAINT  IS  VIOLATED,  CR  SET 
CONSTRAINT    IS    VIOLATED. 


IF    NC     IMPLICIT 

KE»1     IF    ANY     IMFLICI" 


1 1  ME  N  S I  ON 
XI  =  X(  1  ) 
X2  ■  X(2) 
KE  -  0 
X(3 )  =  XI 
IF    (X(3» 


X(4) 


¥    1.732051*X2 

LT.    0.     .OR.     X(3)     .GT.    6.)     GO    TO    1 
X(4)    =    Xl/1. 732051  -X2 
IF    (X (4)     .GE.    0.  )        RETURN 


KE    =    1 

RETURN 
END 


FUNCTION 
DIMENSION 


FE(XJ 
X  (4) 


THIS    IS    THE    OBJECTIVE    FUNCTION. 


FE=    -(  X(2) 

R  ET'JR  N 
END 


**3    * (9.-( X( l)-3. )**2 )  / (46.76  533)  ) 
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C  PLACEMEN"1"    PROCESS     IS    REPEATED    UNTIL    THE    VERTEX    HAS  C 

C  BECOME    FEASIBLE.       IF    THIS    FAILS    TC"  HAPPEN    AFTER  C 

C  5*N+10    DISPLACEMENTS!     THE     SOLUTION    IS    ABANDONED.  C 

C  AP^ER    EACH    VERTEX     IS    ADDED    TO    THE    COMPLEX,     THE  C 

C  CURRENT    CENTROID    IS    CHECKED    FOR     FEASIBILITY.       IF  C 

C  IT    IS    INFEASI3LE,    THE    LAST    TRIAL    VERTEX    IS    ABANO-  C 

C  ONED     AND    AN    EFFORT    TC    GENERATE     AN    ALTERNATIVE  C 

C  TRIAL     VERTEX     IS    MADE.       IF    5*N+1C    VERTICES    ARE  C 

C  ABANDONED    CONSECUTIVELY,    THE    SOLUTION     IS    TER-  C 

C  MINATED.  C 

C  C 

C  IF    AN    INITIAL    COMPLEX    IS    ESTABLISHED,     THE    BASIC  C 

C  COMPUTATION    LOOP    IS    INITIATED.       THESE    INSTRUCTIONS  C 

C  FIND    THE    CURRENT    WORST    VERTEX,     THAT    IS,     THE    VERTEX  C 

C  WITH    THE    LARGEST    CORRESPONDING    VALUE    FOR    THE    OB-  C 

C  JECTIVE    FUNCTION,     AND    REPLACE    THAT    VERTEX    BY    ITS  C 

C  OVER-REFLECTION    THROUGH    THE    CENTROID    OF    ALL    OTHER  C 

C  VERTICES.       (IF    THE    VERTEX    TO    BE    REPLACED    IS    CON-  C 

C  SIDERED    AS    A    VECTOR     IN    N-SPACE,     ITS    OVER-  C 

C  REFLECTION    IS    OPPOSITE    IN    DIRECTION,     INCREASED    IN  C 

C  LENGTH    BY    THE    FACTOR    1.3,     AND    CCLLINEAR    WITH    THE  C 

C  REPLACED    VERTEX    AND    CENTROID    OF    ALL    OTHER  C 

C  VERTICES.)  C 

C  c 

C  WHEN    AN    OVER-REFLECTION     IS    MOT    FEASIBLE    OR    REMAINS 

C  WORST,    IT    IS    CONSIDEREC    NOT-P ER V IS3 I SLE    AND     IS  C 

C  DISPLACED    HALFWAY    TOWARD    the    CENTROID.       AFTER    FOUR  C 

C  SJCH     ATTEMprs     ARE    MACE   UNSUCCESSFULLY,     EVERY    FIFTH  C 

C  ATTEMPT    IS    MADE     BY    REFLECTING    THE    OFFENDING    VERTEX  C 

C  THROUGH    THE   PRESENT    BEST    VERTEX,    INSTEAD    OF  C 

C  THROUGH    THE    CENTROID.       IF    5*N+10    DISPLACEMENTS    AM C  C 

C  OVER-REFLECTIONS    OCCUR    WITHOUT    A    SUCCESSFUL  C 

C  PERMISSIBLE)    RESULT,    THE    CURRENT    BEST    VERTEX    IS  C 

C  TAKEN    AS     AN    INITIAL    FEASIBLE    POINT    FOP    A    RESTART  C 

C  RUN    OF    THE    COMRLETE    PROCESS.        RESTARTING    IS    ALSO  C 

C  UNDERTAKEN   WHEN    6*NV+10    CONSECUTIVE    TRIALS    HAVE  C 

C  BEEN    MADE    WITH    NO    SIGNIFICANT    CHANGE    IN    THE    VALUE  C 

C  OF    THE    OBJECTIVE    FUNCTION.        IN    ALL    CASES,    RESTART-  C 

C  ING    IS    INHIBITED    IF    THE    LAST    RESTART     CIO    NOT    PRO-  C 

C  DUCE    A     SIGNIFICANT     IMPROVEMENT     IN    THE     VINIMiJM  C 

C  ATTAINED.  C 

C  C 

C  IT    IS    RECOMMENDED     THAT    THE    USER    READ    THE    REFER-  C 

C  EMCE    FOR    FURTHER    USEFUL     INFORMATION.     IT    SHOULD    BE  C 

C  NOTED    THAT    THE    ALGORITHM    DEFINED    THERE    HAS    BEEN  C 

C  ALTERED    TO    FIND    THE    CONSTRAINED     MINIMUM,     RATHER  C 

C  THAN    THE     MAXIMUM.  C 

C  C 

C  c 

C  REMARKS  C 

C  C 

C  THE     INTEGER    PROGRAMMING    OPTION    WAS    ADOED    TO    THIS  C 

C  PROGRAM    AS    SUGGESTED     IN    REFERENCE    (2).       A    MIXED  C 

C  INTEGER/CONTINUOUS     VARIABLE    VERSION    CF    8CXPLX  C 

C  WOULD    BE    EASY    TO    CREATE    BY    DECLARING    "MP"    TO    BE    AN  C 

C  ARRAY    OF    NV    CONTROL    VARIABLES    WHERE    IP(I)=1    WOULD  C 

C  INDICATE    THAT    THE    I-TH    VARIABLE    IS    TO    BE    CONFINED  C 

C  TO    INTEGER    VALUES.        EACH    STATEMENT    DF    THE    FORM  C 

C  ' 1=     (I P    .EQ.     1) '     ETC.    WOULD    THE\    NEED    tq    BE    AL-  C 

C  TERED    TO     'IF    < I P ( I )     .EQ.     D'     ETC.,    WHERE    TH E    SUB-  C 

C  SCRIPT    IS     APRRO°R  I  AT  ELY    CHOSEN.       NORMALLY,     XU     AND  C 

C  XL    VALUES    ARE    ALTERED    TO    BE    AN     E^SILON     ♦WITHIN'  C 

C  ACTUAL    VALUES     DECLARED    BY    THE     USER.        THIS    ADJUST-  C 

C  MENT     IS    NOT    MAOc    WHEN     IP»i.  C 

C  C 

C  NOTE:        NO    NON-LINEAR    PROGRAMMING    ALGORITHM    CAN 

C  GUARANTEE    THAT    THE    ANSWER    =OUMD    IS    THE    GLOBAL 

C  MINIMUM,     RATHER     THAN     JUST    A     LOCAL    yINIMLM.       HOW-  C 

C  EVER,     ACCORDING    TO    REF.    2,     THE    CDM°i_EX    "ETHOD    HAS  C 

C  AN    ADVANTAGE    IN    THAT     IT    TENDS    T0    CIND    THE    GLOBAL  C 

C  minimum    MQPE    FREQUENTLY    THAN    MANY    OTHER    NCN- 

C  LINEAR    PROGRAMMING    ALGORITHMS.  C 
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c  ** 

C  IT    SHCULD    BE    NOTED    ThAT    THE    AUXILIARY    VARIABLE  C 

C  FEATURE    CAN    ALSO    3E     USED    TO    DEAL    WITH    PROBLEMS  C 

C  CONTAINING    EQUALITY     CONSTRAINTS.       ANY    EQUAL  TIY  C 

C  CONSTRAINT    IMPLIES    THAT    A    GIVEN    VARIABLE    IS    NOT  C 

C  TRULY     INDEPENDENT.       THEREFORE,     IN    GENERALt    ONE  C 

C  VARIABLE    INVOLVED    IN    AN    EQUALITY    CONSTRAINT    CAN  C 

C  BE      RENUMBERED    FROM    THE    SET    OF    NV    INDEPENDENT  C 

C  VARIABLES     AND    AODED    TO    THE    SET    CF    NAV    AUXILIARY  C 

C  VARIABLES.       THIS    USUALLY    INVOLVES     RENUMBERING  C 

C  THE    INDEPENDENT    VARIABLES    OF    tHE    GIVEN    PROBLEM.  C 

C  C 

C  SLBROUTINES    AND    FUNCTICNS    REQUIREC  C 

C  C 

C  SUBROUTINE    'BOUT'     ANC    FUNCTION     'FBV     ARE    INTEGRAL  C 

C  °ARTS    OF     THE    BOXPLX    PACKAGE.  C 

C  C 

C  TWO    FUNCTIONS    MUST    e5   SUPPLIED    eY    THE    USER.       THE  C 

C  FIRST,    ,KE(X),    IS    USED    TO    EVALUATE    THE     IMPLICIT  C 

C  CONSTRAINTS.       SET    K  E=0    AT    THE    BEGINNING    OF     THE  C 

C  FUNCTION,    THEN    EVALUATE    THE     ILLICIT    CONSTRAINTS.  C 

C  IN    THE    EXAMPLE    ABOVE.    THE    FIRST    CCNSTRAINT,    X(2),  C 

C  MUST     BE    WITHIN    THE    RANGE     (0.     .L  E .    X(3J     .LE.    6).  C 

C  THE     SECONC    CONSTRAINT    X(4),     MUST    BE    .GE.    0.     .       IF  C 

C  EITHER    CONSTRAINT    IS    NOT    WITHIN    THESE    3CUNDS,  C 

C  CCNTROL    IS    TRANSFERRED    TO    STATEM=NT    1,     ANC    KE     IS  C 

C  SET    TO    "1  "    AND    CONTRCL    IS    RETURNED    TO    BOXPLX.  C 

C  C 

C  THE    SECONO    FUNCTION    THE   USER    *UST    PROVIDE    EVALU-  C 

C  ATES    THE    OBJECTIVE    FUNCTION.        IT    IS    CALLEC    FE(X)  C 

C  AS     SHOWN    IN    THE    EXAMPLE    ABOVE,     AND    FE    MUST    BE    SET  C 

C  TO    THE    VALUE    OF    THE    CBJECTIVE    FUNCTION    CQRRES-  C 

C  PONDING    TO    CURRENT    VALUES    OF     THE    NV    INDEPENDENT  C 

C  VARIABLES     IN    ARRAY     *X  •  .  C 

C  C 

C  REFERENCES  C 

C  C 

C  BOX,       M.J.,    A    NEW    METHOD    OF    CONSTRAINED    CPTIMI-  C 

C  ZATION    ANO    A    COMPARISON    WITH    O^HER    METHODS"*  C 

C  COMPUTER    JOURNAL!    3     APR.    '65,     PP.    42-52.  C 

C 

C  BEVERIDGE    G.,    AND    SCHECHTER    R.,    "OPTIMIZATION:  C 

C  THEORY    ANC    PRACTICE",    MCGRAW-HILL.     1970.  C 

C  C 

C  .  C 

C  PROGRAMMER 
C 

C  9.R.     HILLEARY    1/1966. 

C  REVISED    FOR    SYSTEM     260    V1967 

C  CORRECTEO      1/1969 

C  REVI SED/EXTENOED    BY    L  .NO  LA N/R .HI LLEAR Y       2/1975 

C  CORRECTED      3/1976 
C 

c 

IMPLICIT    REAL  *8     (A-H.O-Z) 

REAL*B    V(50  ,50)  ,FUN(  50)  , SUM  (25)  , CEN(25)  ,XS(NV  ),BU(NV), 
*BL(NV) 

C 

KV    =     5 

EP    =    l.D-6 

NTA    =    2000 

IP    (NTZ.GT.O)    NTA    =    NTZ 

R    s    R7 

IF(R.LE.O.0O.OR.R.GE.l .CO )R  =  1 . DO/ 2  .DO 

N  VT    =    N  V+NA  V 
C 
C  TOTAL    VARS,     EXPLICIT    OLUS     IMPLICIT 

NT   =     0 
C  CUR^  ENT    T^ IAL    NO. 

NFT    =    0 
C  CURRENT    NO.     OF     PERMISSIBLE    tqialS 

NTFS    =    0 
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C  CURRENT  NQ •  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 

C 

C  CHECK    FEASIBILITY    GF    START    POINT 


C 


C 

C 


DC   4    1=1, NV 

VT    =     XS(I ) 

IF    (BL(  D.LE.VT  )    GO    TO    1 

II    =    -I 

VT    =    3LI II 

GO   TO    2 

1  IF    (BUm.GE.VT  )    GO   TC    3 
II    =    I 

VT    =    BU(  I) 

2  IF    (NPR.GT.O)     WRITE    (6t49)     II 

3  V(I,1)     =   VT 
CEN(I )    -    VT 

IF    (  IP.EQ.l  )    GO    TO   4 

3L(I  )=BL(I  )+OMAXi(EP,  EP*OABS(BL  (I  )  )  ) 

3U(I)=BU( I )-OMAXl (  cPi  EP*CABS ( BU  (  I  ))  ) 

4  SLM  (I  )     =    VT 


NUMBER    OF    CONSTRAINT    EVALUATIONS 
I     =    1 

IF    (Kc ( V<1,1) ). EQ.O)     GC    TO    5 
I F    (NPR  .LE.O)    GO    TO    12 
WRITE     (6,50) 
GQ    TO    12 

5  NFE   =    1 
C 

C  NUMBER    OP     VERTICES    (  K)     =    2    TIMES    NO.    CF    VARIABLES. 

K    =    2*NV 
r 

C         NUMBER    OF    DISPLACEMENTS     ALLOWED. 

NLIM    =     5*NV>10 
C 

C  NUMBER    CP    CONSECUTIVE    TRIALS    WITH    UNCHANGED    FE   TO 

C        TERMINATE. 

NCT    =    NLIM+NV 

AL°HA=1  .3D0 

FKM=FK-1  .DO 
3ETA=ALPHA>1.D0 

C 

C  INSURE    SEED    OF    RANDOM    NUMBER    GENEcATO'    IS    ODD. 

IGR    =    R*1.0  7 

IF    (MOO<  IQR  ,2)  .S3  .J)     IQR-IQR^lOl 

r 

C  SET    UP    INITIAL    VERTICES 

FUN  (I  )    =    FEtVdfill 
Y"N-  =    FUN(l) 

6  FI    =    l.DC 
FUNOLD    =    FUM1I 

r 

DO    15    I»2»K 

FI  -.    Fin.oo 

LIMT    =    0 

7  LIMT    =    LIMT+1 
C 

C  END    CALCULATION    I?     FEASIBLE    CENTROIC    CANNC"     BE     FCUNC. 

I  F    (L  IMT  .GE  .NLI  1)    GO     TO    11 
C 

DO    8    J=l  ,NV 
C 
C  RAN  COM    NUMBER    GENERATOR       (RANDU) 

I  CR    -    I  GR*65539 

IF  (IQR.LT.O)  I  3R  =  I QR+ 21 47432647+1 

RCX  =  I QR 

RCX  =  RQX*. 465661 30-9 

V(Jrl)  =  BL  (J  )+RQX*(  BU(  J  1-3L(  J  M 

IF(  (V(  J,  I  )+.5DO  )  .  GT.  214  7  43  3647  .)WRI7  =  {6,100) 
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IF    (IP.EQ.l )V(J, I )=IDINT(V( J, I 1+.5D0) 
S    CONTINUE 
C 

OC    10    L-1,NLIM 

NC5    =    NCE+l 

IF    (KE< V<1,  I  )  >.EQ.O)    GO    TO     13 
C 

00    9    J=1,NV 

VT   =  (V(  J,I)KEN{  J)  )*.5D0 

IF(  (  VT  +  .5D0  ).GT.214743  3647.)WRITP(6,  100  ) 

IF    (IP.EQ.l)     VT-IDINT    (VT+.5D0) 

V(Jf  I  )    «    VT 
9    CONTINUE 
C 

10  CONTINUE 
C 

11  IF    (NPR  .LE.O)    GO    TO    12 
WRITE     <6,51)    I 

CALL    BOUT     (NT,NPT  ,NFE,\CE,  NV,  MVT,V,I  ,FUN,CEN,  I) 

12  IER    =    -1 
GC   tq   48 

r 

13  DO    14    J=1,NV 

SUM(J)     =    SUM (J) +V(J, I ) 

14  CEN( J)      -    SUM( J) /FJ 
C 

C         TRY    TO    ASSURE    FEASIBLE    CENTROID    FOR    STARTING, 
NCE    -    NCS+i 

IF    (<E(  CEN).EQ.  0)    GO     TO    60 
SUM(J)     =    SUM(J)    -V(  J,  I) 
GO    TO    7 
60    NFE   =    NFE+1 

FUN<  I)     -    FE(V(1,  I  )  ) 

15  CONTINUE 
C 

C  END    CF    LOOP   SETTING   OF     INITIAL    COMPLEX. 

IF    (  NPR.LE.  0)    GO     TO    17 

CALL     30UT     (NTtNPT  ,NFE,  NCE,NV,M  VT  ♦  V ,  K,FUNTCEN,  0) 
C 
C  PINO    THE     WORST    VERTEX,    THE    'J'TH. 

J    =    1 
C 

HO    16    1=2, K 

IF    (FUN (J ) .SE.FUN(  I  )  )    GC    T0    16 

J    =    I 

16  CONTINUE 
C 

C  BASIC    LOOP.    ELIMINATE    EACH    WORST    VERTEX    IN    TURN.     IT 

C         MUST   BECOME    NO    LONGER    WORST, NOT    MERELY    I-MPRCVED.    FIN 
C  N£XT-T0-WORST    VERTEX,  THE     »JN,TH    ONE. 

17  JN»1 

IF    (J.EQ.l)     JN    =    2 
C 

DC    IB     I  =1,K 

IF    (  I  .EQ.J  )    GO    TO     13 

IF    (FUN  (JN  )  .GE.FUN(  I  )  )    GO    rj     iq 

JN   '    I 
13    CONTINUE 
C 

C  LIMT=NUMBER     OF    MOVES    DURING    THIS    TR  I  &L    TOW  ARC    TrE 

C         CENTROID    DUE     TO    FUNCTION    VALUE. 

LIMT    =    1 

C         COMPUTE    CENTROID    AND    OVER    REFLECT    WCRS^    VERTEX. 
C 

DO    19    I  =1  ,NV 

VT    =     V( I ,J) 

SUM(  I  )     a    SUMU  J-VT 

CEN(  I  )      =     SUM  (I)  /f=KM 

VT    =    5ETA*CEN(  I  )-ALPHA*^T 

IF(  (VTf.500  ).*T  .2  147  43  3  647.  )w  R  ITE  (  6  ,  1  JO  ) 

IF    (IP.EQ.l)     VT    «rOINT<  VT4-.5D0  ) 
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c 
c 


c 

c 
c 


c 
c 
c 


c 
c 
c 


r 

f 


c 
c 
c 
c 
c 
c 


INSURE    THE    EXPLICIT    CONSTRAINTS     APE    OBSERVED. 

19  V(I,J)     =    DMAXKOMIN1  (VT,BU(I)  )  ,BL(I)  ) 

NT   =    NT+1 
CHECK    FOR    IMPLICIT    CONSTRAINT    VIOLATION. 

20  OC    25     N=1,NLIM 
NCE    =    NCE+i 

IF    C<E(V(1,  J)  )  .EQ  .0)    GO    TO    26 

EVERY    'KV'TH    TI  ME , OVER-REFLECT    THE    CFFENDING    VERTEX 
THROUGH    THE    BEST    VERTEX. 

IF     (MOD(N,KV)  .NE.OI    GO    TC   22 

CALL    FBV    (K, FUN, Ml 

DC  21     I»lfNV 

VT    =    BETA*V(I,M)-ALPHA*V(I,J)  ' 

IF<  (VT  +  .5D0  ).GT .2147432647. )WR I TE(  6,100) 

IF    (IP. EC!)    VT    »ID1NT(VT+.5D0) 

21  V(I,J)    =   DMAXKDMINKVT, 8U(I)  ), BLUM 


GC   TO   24 
CONSTRAINT    VIOLATION: 


MOVE    VJEW     POINT    TOWARD    CENTROID. 


22  00  23    1*1. NV 

VT    s  (CEN(I)-H/  (I,  J  )  )*.5D0 

I=( (VT+.5D0).GT.2147  46  3  647. ) WR IT 5 (6 ,100 ) 

IF    (ID. ECU     VT    -IDINT(  VT*.3DQ) 
V  (I, J)     =*    VT 

23  CONTINUE 

24  NT    =    NT+i 

25  CONTINUE 

I  ER   =    1 

CANNOT    GET    FEASIBLE    VERTEX    3Y    MOVING    TOWARD    CENTRCID, 
OR    BY    OVER-REFLECTING    TH<?U    THE    3  EST    VERTEX. 

IF    (MPR.LE.OI     GO    TO    42 

WRITE    (6,52)    NT,  J 

CALL    BOUT     (NT,NP7 ,NFE,MCE, NV,MVT, V, K, FUN ,C5M, J) 

GC   TO    42. 

FEASI3LE    VE*TEX    FOUND,     EVALUATE   THE    03JECTIVE    cUNCTION. 

26  NFE    =    NFE+L 
FUMTRY    -    FE( V( 1, J ) ) 

TEST    TO     SEE    IF    FUNCTION    VALUE    HAS     NCT    CHANGED. 
AFO    =DABS(  PJNTR/-FUNOLD  ) 
APXOMAX!  (DABS  ( £P*FUNCLC) ,  EP) 

ACTIVATE    THE    FOLLOWING    TWO    STATEMENTS    FOR    DIAGNOSTIC 
PURPOSES    ONLY. 

WRITE (6, 99) J, AFO, AM Xf FUNTR Y , FUNGLC , FUN < J ) ,FLN ( JN) , 
*NTFS ,N 
99    FCRMAT     (IX, 13, 6515. 7, 215) 

IF    ( AFO.  GT  .amx)    GO    TO     27 

NTFS    =    NTFS+l 

IF    (N7FS.LT. NCT)     GO    TO    23 

I  ER    =    0 

IF    (NPR.LE.O)     GO    TO    42 

WRITE     (6,53)     K 

CALL    BOUT     (NTfNPT,NFE,NCE,NV,NVT, V,K,FUN,CEN, C) 

GO    TO     42 

27  NTFS    =     0 

IS    THE    NEW    VERTEX    NO    LCNGE'    ^ORST? 

28  IF    (FUNTRY.LT.FUN (JN) )     GO    TO    34 


74 


C         TRIAL    VERTEX     IS    STILL    WORST;    ADJUST    TOWARD    CENTRGID. 
C  EVERY    'KV'TH    TIME,    QVER-R£CLECT    THE    C?CENQING    V  ERT  EX 

C         THROUGH    THE    BEST    VERTEX. 

LIMT    =    LIM7+1 

IF    (MOO (LIMT, KV) .NE.O)     GO    TQ    30 

CALL     FBV    (K.  ,FUNfM  ) 
C 

DO    29    I =  1,NV 

VT    =    BETA*V(I,M)-ALPHA*V(If J) 

IF(  (VT  +  .5CO  ).  GT.  2147483  647  •  )W  RITE  (6,100) 

IF    (IP.EQ.l)     VT    aXOINTt  VT+.530  ) 
29    V(I,J)     =    DMAXKDMINH  VT,3U(I  )  )  ,3L(I  )) 


GO   TD    32 
C 

20  00    31    I-1,NV 

VT   =  (CEM(  I)+V(  I,  J  )  M.5D0 

IF( ( VT +. 3 00 ).GT •214748  36 47. >WR IT  6(6, 100  ) 

I11    (IP.EQ.l)     VT    *IDINT{VT+.5DG) 

V(  I, J  )     =    VT 

21  CONTINUE 

r 

32  IF    (LIMT.LT.NLIM)    GO    TO    33 

C  CANNOT    MAKE    THE    •  J«  TH    VERTEX    NO    LCNGE3    WORST     3Y 

C         DISPLACING   TOWARD    THE    CSNTROID    OR     BY    OVER-REFLECTING 

C  THROUGH    THE    BEST    VERTEX. 

IER    =    2 

IF    (NPR     .LE.    0)       GO    TO    42 

WRITE     (6,5  2)        NT,     J 

C  *LL    BOUT    (  NT,NPT,NFE  ,NCE,N V,N VT  ,  V  ,  K,FUN  ,CEN,  J) 

GO   TO    42 

33  NT    =    NT+i 
GO    TO    20 

C 

C  SUCCESS:       WE    HAVE    A    REPLACEMENT    FOR    VERTEX    J. 

34  FUN(J)  *  FUNTRY 
FUNOLD  =  FUNTRY 
NPT    =    NPT+1 

C 

C         EVERY   lOO'TH    PERMISSIBLE    TRI AL»R ECQVPUTE    CENTROID 

C  SUMMATION    TO     AVOID    CREEPING    ERROR. 

IF     (MHDtNPT, 100  )  .NE.O  )    GO    TO     37 
C 

c 
c 
c 
c 


c 
c 

c 

r 


DC   3  6    I=l,NV 
SCM( I )     =    O.DO 

DO    35    N«1.K 
25    SJM(  I  )    =    SUM(  I)*V( I,N) 

CEN(  I  )     =    SUM(I)  /FK 

36    CONTINUE 

LC    =    0 
GO    TO    39 

27    DO    38    I  =l,NV 

38    SUM(I)     =    SUM(  I)+V(  I, J) 

LC    =    J 

3°     IF    (N^R  .LE.O  )    GO    TO    40 

IF    (MOD ( NPT tNPR) . NE.O)     GC    T0    40 

CALL     BOUT     (NT  ,NPT  ,  NFE  ,  NCE,  N  V,  N  VT  ,  V  ,  K  ,F'J  N  ,C  EN,  LC  ) 


C  HAS    THE    MAXIMUM    NUMBER     CF     TRIALS     BEEN    REACHED    WITHCU' 

C  CONVERGENCE?     IF    N1T,     GO    TO    NEW    TQIAL. 

40    IF    (NT.GE.NTA)    GO    TO    41 

C  NEXT-TG-WO£;ST    VERTEX    NOW    BECOMES    WORST. 

J    =    JN 
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GO   TO    17 

41  IER    =    3 

IF    (NPR  .GT  .0)    WR  ITE    (  6t54) 
C 

C         COLLECTOR    POINT  FOR    ALL    ENDINGS. 

C  1)       CANNOT    DEVELOP    FEASIBLE    VERTEX.  IER    =    1 

C  2)       CANNOT    DEVELOP    A    NC-LCNGEP-WORST    VERT  EX  •       I  "R    =    2 

C  3)       FUNCTION    VALUE     UNCHANGED    FOR    K    TRIALS.  IER    =    0 

C  4)       LIMIT    ON    TRIALS    REACHED.  IER    =     3 

C  5)       CANNOT    PINO    FEASIBLE    VERTEX    AT     START.  IER    -    -1 

42  CONTINUE 
C 

C         FIND   3  EST    VERTEX. 

C^LL     FBV    (K,FJN,M) 

IF    (  IER.GE.3  )    GO    TO    44 
C 

C         RESTART     IF    THIS    SOLUTION    IS    SIGNIFICANTLY    BETTER    THAN 
C         THE     PREVIOUS, OR    IF    THIS     IS    THE    FIRST    TRY. 

IF    ( NPR.LE.O)    GO     TO    43 

WRITE    (6,55)    (MtYMNfFUN(.M)  ) 

43  IF     (FUN(M> .GE.YMN)    GO    TC    47 
IF(DA3S(FUN(M)-YMN)  .LE.DNAXi(EP,EF*YMNn     GO    TO    47 

C         GIVE    IT    ANOTHER    TRY    UNLESS    LIMIT    ON    TRIALS    REAChEC. 

44  Y  *N    =    F  UN  i  M  ) 
FUN(l)     =     FJN(M) 

C 

DO    45    I =1,NV 
C6N(  I  )     =    V(  I,M) 
SUM (I )     =    V (  I,M) 

45  Vt 1,11     =    V(  I ,M) 

DC    ^6     1=1,  NVT 

46  XS(I  )     =    V(  I  ,M) 

IF     (IER.LT.3)     GO    TO    6 

47  IP    ( MPR.LE. 0)    GO     TO    43 
CALL     BOUT     (NT, NOT  ,NFE  ,  NCE,  NV,  N  VT  ,  V  ,K  ,FUN  ,  V  (  1,M)  ,-D 
WRITE     (6,56)     FUN(M) 

48  RETURN 

49  FORMAT  ( 'OlMuEX     AN  C    DIRECTION    OF    OUTLYING    VARIABLE1, 

*  •     AT    START'   ,1  5) 

50  FORMAT  {  »0  IMPLICIT    CONSTRAINT    VIOLATED    AT    START.  t? 
*OX,'OEAD    END.») 

51  FOPMAT(  'OCANMOT    FIND     FE A  SI BLE '  , 14  , 
*«TH    VERTEX    OP    CENTROIC    AT    START.1  J 

52  =ORMAT('GAT  TRIAL  ',14,'  CANNOT  PINC  FEASIBLE  VERTEX1, 
*' WHICH  IS  NO  LONGER  WOR ST'  , I  4, 1 5X  , 'RE  STAR T  FROM  BEST', 
*'  VERTEX.'  ) 

53  =OP^AT  (40HOFUNCTION  HAS  3EEN  AL^CS"  UNChdNGEC  FOR  15, 
*7H    TR  IALS) 

54  FCRMAT     (27H0LIMIT     CM    TRIALS    EXCEECEl.     ) 

55  FORMA  T(  •  08  EST    VERTEX    IS    NO  .  •  ,  I  3  ,  »        OLD     *IN    WAS     ', 
♦E15.7,  •       NEW    MIN    IS     '  ,E15.7) 

56  =GRMAT    COMIN    03JECTIVE    PJNC^ION     IS     SE15.7) 
100    FORMAT (  '0'  , 'TRUNCATION    ERROR    IN     BCXPLX') 

END 

SUBROUTINE    FBV     (K,FUN,*) 
IMPLICIT    REAL*3  (  A-H,0-Z) 
P.EAL*8    FUN(  30) 
M    =    1 

-)  n    i     1=2    K 

IF     (FUN(  M)  .LE.FUN  (  I  )  )     GO    T3     \ 
M    =    I 
1    CONTINUE 
RE^U'N 
END 


SUBROUTINE    30UT    ( N T, NP T ,NFE ,NC 5 , N \ ,NVT ,V ,K ,FN ,C , I K ) 
IMPLICIT    REAL*8    (A-HfO-ZJ 
REAL*8    V<  50,50)  ,FN(50)  ,C(25) 
WRITE     (6,4)    NT,NPT,NFE,NCE 

DC    1    1=1, K 

WRITE     (6,5)     FN(  I)  ,  (  V(  J, I  )  ,J=1  ,NV) 

IF    (NVT.LE.NV)     GO    TO    1 

NVP    =    NV+1 

WRITE    (6,6)     (V(  J,  I),  J=NVP,NVT) 

CONTINUE 


IF    {  IK.NE.O  )    GO    TO    2 
C 

WRITE    (6,7)     (C(I)  ,1=1  ,NV) 

RETURN 

2  IF    (  IK.GE.O  )    GO    TO    3 
WRITE     (  6,3)     (C(  I)  ,1=1  ,NV) 
RETURN 

3  WRITE     (6,9)     IK*  (C(I)t  I*LfNV1 
RETURN 

C 

4  FCRMATCONO.    TOTAL    TRIALS    =    ',I5,4X, 
*'N0.     FEASIBLE    TRIALS     =    ',I5,4X, 
*»N0.    FUNCTION    EVALUATIONS    =    '  ,  I  5,  4X, 
*'N0.    CONSTRAINT    EVALUATIONS    =     ',15/, 
*«Q  FUNCTION    VALUE', 6X, 

*« INDEPENDENT    7  ARI  ASL  ES  /  C  EP  ENDENT    OR     IMPLICIT', 
*'CCNSTRAINTSf  ) 

5  FORMAT     (  ih    ,E13.7,2X, 7E14. 7/( 21X,7E14.7  )  ) 

6  FCRMAT     (2iX,7E14.7) 

7  FORMAT    (10H0CENTRQID    11X ,7E14. 7/ ( 21X , 7E14. 7  )  ) 

8  FOR*AT     CO       BEST    VER  TE  X«  f  7X ,  7E  1  A.  7/  (  2IX  ,  7E  1 4.  7)  ) 

9  FCRMAT     COCENTRQID    LESS    VX • , 12 , 2X , 7E14. 7/ ( 21X, 7E14.7) ) 
END 

C  C 

C  A    /*    CARC    GOES    HERE  C 

C  C 

//GO.SYSIN    DO    * 

i  2  ?  1  2  2  110 

0.0  0.005  15.  C~ 

13  50  2 

100.0  0.4  10.0  1.0  0.0  STEP 

14  50  2 

16.0  0.6  4.0  1.0  0.0  STEP 

18 
13    13       1100.0  3.0  7.63636262 

2  3    14      -400.0  3  .0  4.  3 

3  2    2    5      -32.0  4.3  7.0 

4  2    2    6      24.  C  4.3  5.73333323 

5  117      100  .0 

6  1    2     3      16.0 

7  2  3  7  1.0  0.0 
3  2  4  3  1.0  0.0 
9    2    5    7      1.0                    0.0 

10  2    6    S      1.  0  0.0 

11  2    7    9      1.0  7.0 

12  2    31  0      1  .0  4.3 

12  2    911      1.  0  12.0 

14  21012      1  .0  2.0 

15  112    9      2.0 

16  11 110      4.0 

17  lii    1      -1.0 

13  112    2      -1.0 

1.0  0.0  S  TE  P    1 

1,0  0.0  STEP    2 

MOW  RE Y    10    JAN    197  3 
THESIS    SIMULATION       Yl 
"OWREY     10    J  AN     1973 
THESIS    SIMULATION       Y2 
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